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APPLICATION  OF  HIERARCHICAL  HIGHER-ORDER  TANGENTIAL  VECTOR 
FINITE  ELEMENTS  IN  A  HYBRID  FEM/MOM  METHOD 

Hao  Wang,  ChunLei  Guo,  and  Todd  H.  Hubing 
Department  of  Electrical  and  Computer  Engineering 
University  of  Missouri-Rolla 
Rolla,  MO  65409 


ABSTRACT 

Hybrid  FEM/MoM  methods  combine  the  finite 
element  method  (FEM)  and  the  method  of  moments 
(MoM)  to  model  inhomogeneous  unbounded 
problems.  These  two  methods  are  coupled  by 
enforcing  field  continuity  on  the  boundary  that 
separates  the  FEM  and  MoM  regions.  Hierarchical 
higher-order  tangential  vector  finite  elements 
(TVFE’s)  are  of  practical  interest  because  they  can  be 
easily  combined  with  low-order  elements  to  improve 
the  accuracy  of  numerical  solutions.  This  paper 
presents  a  hybrid  FEM/MoM  formulation  applying  a 
set  of  hierarchical  TVFE’s  developed  by  Webb  and 
Forghani.  Higher-order  FEM  elements  are  coupled  to 
MoM  elements  based  on  Rao- Wilton-Glisson  (RWG) 
functions.  The  FEM  matrix  assembly  procedure  is 
described  in  sufficient  detail  to  aid  other  investigators 
who  wish  to  develop  codes  employing  this  technique. 
Three  practical  electromagnetic  problems  are 
presented  that  demonstrate  the  advantages  of  the 
higher-order  elements. 

I.  INTRODUCTION 

The  hybrid  finite-element-method  and  method-of- 
moments  (FEM/MoM)  can  be  used  to  analyze  many 
kinds  of  electromagnetic  problems  effectively  by 
applying  FEM  to  model  the  fields  in  regions  with 
geometric  complexity  and  using  MoM  to  model 
larger,  simpler  structures  outside  this  region  and  to 
provide  an  accurate  radiation  boundary  condition 
(RBC)  to  terminate  the  FEM  mesh.  Both  the  MoM 
and  FEM  are  powerful  methods,  but  each  of  these 
methods  has  its  own  advantages  and  disadvantages. 
MoM  handles  unbounded  problems  very  effectively 
but  is  less  efficient  when  complex  inhomogeneities 
are  present.  Inhomogeneities  are  easily  handled  by 
the  FEM,  which  requires  less  computer  time  and 
storage.  However,  the  FEM  is  most  suitable  for 
bounded  problems.  Hence,  hybrid  FEM/MoM 
methods  that  combine  MoM  and  FEM  are 
advantageous  for  treating  electromagnetic  problems 
involving  unbounded,  complex  structures. 

Conventional  hybrid  FEM/MoM  codes  employ  linear 
tangential  vector  finite  elements  (TVFE’s).  These 
elements  are  commonly  referred  to  as  Whitney 


elements  defined  by  Nedelec  [1].  Because  the 
functions  do  not  impose  normal  component 
continuity  between  tetrahedra,  they  do  not  produce 
the  spurious  modes  that  can  be  generated  by  using 
node-based  elements  [2].  However,  these  elements 
limit  the  accuracy  of  the  finite  element  solution  since 
they  only  support  a  constant  tangential  value  along 
element  edges  and  a  linear  field  variation  inside  the 
element  (CT/LN).  Thus,  when  electric  fields  in  a 
certain  region  vary  quickly,  the  number  of  tetrahedra 
has  to  be  relatively  high  to  obtain  reasonable 
accuracy.  Higher-order  elements  that  support  non¬ 
linear  field  variations  can  be  used  to  model  rapidly 
vaiying  fields  using  fewer  elements.  One  set  of 
higher-order  basis  functions  for  tetrahedra  supports  a 
linear  tangential,  quadratic  normal  (LT/QN) 
representation  of  the  fields.  Basis  functions  of  the 
next  higher  order  have  a  quadratic  tangential,  cubic 
normal  (QT/CuN)  representation  for  the  fields.  A  set 
of  TVFE’s  is  referred  to  as  interpolator  if  values 
within  the  element  can  be  interpolated  from  node  or 
edge  values.  It  is  referred  to  as  hierarchical  if  the 
lower-order  basis  functions  are  a  subset  of  the  higher 
order  basis  functions.  Webb  and  Forghani  [3],  Savage 
and  Peterson  [4],  Graglia  et  al  [5],  and  Andersen  and 
Volakis  [6]  have  employed  LT/QN  basis  functions. 
The  TVFE’s  presented  in  [4]  and  [5]  are  interpolator 
while  those  presented  in  [3]  and  [6]  form  a 
hierarchical  set  with  the  Whitney  TVFE.  Hierarchical 
sets  of  TVFE’s  allow  for  selective  field  expansion 
using  different  order  elements  in  different  regions  of 
the  computational  domain.  Hence,  for  the  regions 
where  the  fields  var  slowly,  the  lowest  order 
TVFE’s  can  be  employed,  while  for  the  regions 
where  the  fields  var  rapidly,  higher-order  TVFE’s 
can  be  employed.  This  can  save  memor  and  CPU 
time  without  compromising  computational  accuracy. 
Andersen  developed  and  applied  mixed-order 
hierarchical  TVFE’s  in  [7]. 

For  MoM  techniques  based  on  EFIE  formulations, 
Nedelec  [1]  presented  a  general  family  of  divergence- 
conforming  functions  whose  lowest-order  member 
was  a  set  of  CN/LT  basis  functions,  known  as  Rao - 
Wilton-Glisson  (RWG)  or  triangular  rooftop 
functions.  These  functions  are  widely  used  for 
representing  surface  currents  in  EFIE  formulations. 
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In  a  hybrid  FEM/MoM  technique,  hierarchical  sets  of 
TVFE's  are  readily  coupled  to  linear  MoM  boundary 
elements  because  the  coefficients  corresponding  to 
any  higher-order  terms  on  the  boundary  can  be  set  to 
zero  to  enforce  the  continuity  of  the  tangential  fields. 

The  coefficient  matrices  generated  by  FEM/MoM 
codes  may  have  large  condition  numbers.  When 
LT/QN  basis  functions  are  used  in  the  FEM,  the 
condition  numbers  of  the  hybrid  matrix  generally 
become  much  larger.  Savage  [8]  showed  that 
interpolator  vector  basis  functions  are  generally 
better  conditioned  than  hierarchical  vector  basis 
functions.  However,  only  the  condition  numbers  of 
individual  element  matrices  were  studied  in  [8]. 
Andersen  [9]  examined  the  inter-relationships 
between  the  condition  numbers  of  element  and  global 
matrices  based  on  various  interpolator  and 
hierarchical  TVFE’s  using  a  cavity  resonator 
example.  However  because  they  were  solving  for  the 
eigenvalues  of  a  cavity  resonator,  the  condition 
numbers  of  the  global  FEM  matrices  were  not 
considered. 

In  this  paper,  the  hybrid  FEM/MoM  formulation 
using  the  LT/QN  TVFE’s  described  by  Webb  and 
Forghani  in  [3]  is  developed  and  applied  to  different 
electromagnetic  problems.  Section  II  presents  the 
hybrid  FEM/MoM  formulation.  Section  III  presents 
the  FEM  matrix  assembly  procedure  for  LT/QN 
TVFE’s.  Section  IV  presents  a  set  of  numerical 
results  that  demonstrates  the  improved  performance 
of  the  higher-order  TVFE  in  the  context  of  the  3-D 
hybrid  FEM/MoM. 


II.  FORMULATION 

In  the  hybrid  FEM/MoM,  an  electromagnetic 
problem  is  divided  into  an  interior  equivalent  part 
and  an  exterior  equivalent  part.  The  interior  part  is 
modeled  using  the  FEM  and  the  exterior  part  is 
modeled  using  a  surface  integral  equation  method-of- 
moments  technique  (MoM).  The  two  equivalent  parts 
are  coupled  by  enforcing  die  continuity  of  tangential 
fields  on  the  FEM  and  MoM  boundary  [10]. 

2.1  The  Finite  Element  Method  Using  Higher  - 
Order  TVFE’s 


FEM  can  be  used  to  analyze  the  interior  equivalent 
part  by  solving  the  weak  form  of  the  vector  wave 
equation  as  follows  [7]: 


J 

i- 


VxE(r) 

J°>  AoA,  J 


(V 


*  w(r))  +  j  a  ff0£vE(r)  •  w(r) 


dV 


-  J(/7xH(r))» w(r)dS-  Jj’n'(r)»  w(r)*/J’ 
s  r 


where  5  is  the  surface  enclosing  volume  V,  w(r)  is  the 
weighting  function,  and  J'"'  is  an  impressed  source. 
Equation  (1)  shows  that  efficient  finite-element 
analysis  of  electromagnetic  fields  in  3-D  regions 
requires  computation  of  two  element  matrices.  These 
two  matrices  are 

£,,  =  |  Vx  w,  •  Vx wy  dV  (2) 

Fv  =  |,w,  •  Wj  dl’  (3) 

where  w,  represents  the  \th  vector  basis  function  and 

V  indicates  integration  over  one  tetrahedron.  The  six 
edges  and  four  faces  of  a  tetrahedron  are  numbered  as 
indicated  in  Table  1  and  Figure  1  [4], 


Table  1.  Node  and  edge  numbering  scheme  of  a 
_  tetrahedron 


Edge# 

Node  1 

Node  2 

1 

I 

2 

2 

1 

3 

3 

1 

4 

4 

2 

3 

5 

2 

4 

6 

3 

4 

Face# 

Node  1 

Node  2 

Node  3 

I 

1 

2 

3 

2 

1 

2 

4 

3 

1 

3 

4 

4 

2  I 

3 

4 

Figure  1 .  Edge  and  face  definition  of  a  tetrahedron. 

The  linear-tangential,  quadratic-normal  (LT/QN) 
basis  functions  developed  by  Webb  exist  in  two 
forms.  One  is  edge-based  functions,  which  are 
associated  with  tetrahedron  edges.  The  other  is  face- 
based  functions,  which  are  associated  with 
tetrahedron  faces.  The  two  edge-based  LT/QN  basis 
functions  associated  with  edge  i  are, 

w?  =h(Lt]VLn-Lt2VLt))  k  =  1,- - -,6  (4) 

w?  =  ULtlVLtI  +  LaVLn)  k  =  1, -,6  (5) 

where  “el”  represents  the  first  type  of  edge  basis 
function,  and  “e2”  represents  the  second  type  of  edge 
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basis  function.  L  is  the  area  coordinate  associated 

i 

with  the  node  i.  It  is  unity  at  node  i  and  decays  in  a 
linear  fashion  to  zero  at  the  other  three  nodes  of  the 
cell.  it  is  the  length  of  edge  i. 

Figure  2  shows  vector  plots  of  the  edge-based 
functions  in  a  face  of  a  tetrahedron.  The  two  face- 
based  elements  associated  with  face  i  are, 

w f=I,1I,2VLl3-I(1L,3VA2  /  =  1,-4  (6) 

wf2  =  LaLl2VL,3  -  L,2Ll3VLn  /  =  !,•••  4  (7) 

where  “fl”  represents  the  first  type  of  face  basis 
function,  and  “f2”  represents  the  second  type  of  face 
basis  function.  Figure  3  shows  vector  plots  of  the 
face-based  functions  on  a  face  of  a  tetrahedron.  It 
shows  that  the  field  distributions  in  fl  elements  and 
f2  elements  are  similar  but  they  rotate  in  different 
directions. 


^  v  v  \ 

~  ^  v  V  \  V 

*  ^  N  \  \  \  \ 

.  *  %  \  \  \  \  A 

_ . .  ■  1 1 1 1 1 1 


(a) 


(b) 


Figure  2.  Plot  of  the  edge  based  basis  functions. 

(a)  w'1  >  (b)  w'2. 

Using  these  basis  functions,  the  electric  field  E  in  the 
interior  region  can  be  expanded  as  the  sum  of  four 
terms, 

E(r)  =X(^c,wI5  +  Ef  wf  +  £f  wf  +  £fwf  )* 

k=l 

The  basis  function  Wjfc  has  the  following  properties: 

i  (9) 

lk 

et»wf=/;^Y^  =  4-,-4,  (10) 

h 


where  e  is  a  unit  edge  vector  corresponding  to  the 

k1h  edge.  Hence,  the  terms  associated  with  “el” 
elements  can  be  viewed  as  the  main  terms  that 
describe  fields  along  tetrahedron  edges  roughly, 
while  the  terms  associated  with  “ e2 ”  elements  can  be 
viewed  as  adjustment  terms  that  describe  the  field’s 
linear  variation  along  tetrahedron  edges. 


Figure  3.  Plot  of  the  face  based  basis  functions, 
(a)wf,(b)  wf. 


Since  CT/LN  functions  have  one  unknown  per  edge, 
they  generate  6x6  local  matrices.  LT/QN  functions 
have  two  unknowns  per  edge  and  two  unknowns  per 
face  so  they  generate  20  x20  local  matrices. 
Applying  the  LT/QN  basis  functions  to  (1),  a  global 
FEM  matrix  can  be  constructed  as  follows, 


“  jclci 
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jclcl 

Adl 
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The  unknown  coefficients  [En]  are  partitioned  into 
four  types  according  to  their  corresponding  basis 
functions  and  edge  functions.  The  four  categories  are 
interior  edges  of  “el”  type,  which  are  denoted  by  the 
subscript  /,  dielectric  boundary  edges  of  type  “el”, 
which  are  denoted  by  the  subscript  d,  interior  edges 
of  type  “e2”,  and  interior  faces  of  type  “f,  which  are 
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also  denoted  by  the  subscript  ».  The  type  "fl"  and 
type  uJ2n  basis  functions  have  been  combined  into 
one  common  “f  type  because  they  are  essentially  the 
same  when  they  share  the  same  face  and  rotate  along 
the  same  edge.  Ef ,  El'  and  £*  are  set  equal  to  zero 

on  the  MoM  boundary  to  enforce  the  continuity  of 
the  tangential  electric  fields.  Using  this  approach,  the 
MoM  part  (employing  linear  basis  functions)  does 
not  have  to  be  modified  to  work  with  FEM  elements 
of  different  order.  [./*]  is  a  set  of  unknown  complex 
scalar  coefficients  for  the  surface  electric  current 
densities  on  the  FEM  and  MoM  boundary  S.  [g""]  is 
the  source  term,  representing  sources  located  within 
the  FEM  region.  The  elements  of  [A],  [B^],  and  [ginl] 
are  given  by, 


_  r,(Vxw,(r))»(Yxw„(r)) 

7  J* 

V 

+  7<Deoerw.(r),w„(r  )]</(' 

(12) 

=  Jf.(r)»w„(r)rf5 

(13) 

'  + - - - (VxM'") 

J&MoMr 

•  w„(r  )dV 

.(14) 

2.2  The  MoM  Using  EFIE 

The  exterior  equivalent  part  can  be  analyzed  using 
the  EFIE  [8], 

^=Etac(r)+  J[-M(r')xV'G.(r,r') 

L  s 

-  j  ko  r\,  J(r')G.(r,r') + j  V'.  J(r')  V'G,  (r,  r'  )]</S  ( 1 5) 
ko 

The  equivalent  surface  electric  current  J(r)  and 
magnetic  current  M(r)  in  (15)  can  be  discretized 
using  the  Rao-Wilton-Glisson  basis  function  f(r)  [9], 

JM-ftUf-M  (16> 

■-1 

M(r)  =  £(E*)[f.(r)  (,7) 

■-i 

where  Ns  is  the  total  number  of  edges  on  the  FEM 
and  MoM  boundary  S,  and  N d  is  the  total  number  of 
edges  on  the  dielectric  boundary  S d.  E(r)  in  Equation 
(15)  can  be  expanded  using  the  tetrahedral  CT/LN 
basis  function  w''(r)  as  follows, 

E(') =£;(#).  w*>)  0g) 

■-I 

On  the  surface  S,  the  triangular  basis  function  f(r) 
and  the  CT/LN  basis  function  we,(r)  are  related  by, 


wcl(r)  =  nxf(r)  •  (19) 

After  multiplying  by  weighting  functions  f„(r), 
n=l,  ...  N,  the  EFIE  in  Equation  (15)  can  be 
discretized  as  follows, 

[c]k]=[o][r;]-[r].  (20> 

2.3  The  Hybridization  of  FEM  and  MoM 

Equations  (11)  and  (20)  form  a  coupled  and 
determined  system.  Three  different  formulations,  the 
combined  formulation,  the  inward-looking 
formulation  and  the  outward-looking  formulation , 
can  be  used  to  solve  the  coupled  system  [11],  [14], 
The  outward-looking  formulation  was  used  for  the 
examples  in  this  paper.  From  (20), 

JS=C-'DE'J-C-'F'.  (21) 

Substituting  Equation  (21)  into  Equation  (1 1)  yields  a 
determined  matrix  equation, 
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rr 
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A*' 
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jc2c2 

Af 
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4f 

a: 

Ej 

S,n, 

&  ml 


+  [b«C'F‘] 


(22) 


Equation  (22)  can  be  solved  using  iterative  solvers. 
The  preconditioning  technique  reported  in  [14]  can 
be  used  to  improve  the  convergence  rate  and 
accuracy  of  the  iterative  solvers. 


III.  ASSEMBLY 

The  aim  of  the  assembly  procedure  in  FEM  is  to 
construct  the  global  matrix  (11)  by  summing  the 
element  matrix  terms  for  each  tetrahedron  in  the 
mesh  while  guaranteeing  continuity  of  the  tangential 
electric  field  on  the  boundary  between  any  two 
tetrahedra.  For  CT/LN  basis  functions,  the  assembly 
procedure  is  relatively  straight-forward.  However,  for 
LT/QN  basis  functions,  more  details  have  to  be 
considered  in  order  to  get  the  correct  global  matrix. 
This  section  describes  the  assembly  procedure  for 
LT/QN  TVFE’s. 

For  “e/”  elements,  Equation  (4)  and  Figure  2(a) 
indicate  that  the  complex  scalar  £»'  is  the  projection 

of  the  electric  field  onto  the  k,h  edge.  When  the  local 
edge  vector  (as  defined  in  Table  1)  is  reversed, 
Equation  (4)  will  be  reversed  at  the  same  time. 
Therefore,  to  ensure  the  continuity  of  the  tangential 
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electric  field  across  all  edges,  a  unique  global  edge 
direction  must  be  defined  (e.g.  always  pointing  from 
the  smaller  node  number  to  the  larger  node  number). 
Equation  (4)  must  be  multiplied  by  (-1)  if  the  local 
edge  vector  does  not  have  the  same  direction  as  the 
global  edge  direction. 

For  “e2”  elements,  the  continuity  of  the  tangential 
electric  field  across  all  edges  must  also  be  satisfied. 
From  Figure  2  (b)  and  Equation  (5),  it  is  clear  that 
when  the  local  edge  vector  is  reversed,  there  is  no 
change  in  Equation  (5).  Therefore  Equation  (5) 
should  not  be  multiplied  by  (-1)  even  when  the  local 
edge  vector  does  not  have  the  same  direction  as  the 
global  edge  direction.  When  a  FEM  edge  is  on  the 
boundary  between  FEM  and  MoM,  Ef  associated 

with  this  edge  is  set  to  zero. 

For  77”  and  “£2”  elements,  the  continuity  of  the 
tangential  field  needs  to  be  enforced  across  all  faces. 
From  Figure  3,  two  local  e[  can  be  regarded  as  a 

common  global  unknown  only  if  they  share  the  same 
face  and  rotate  along  the  same  edge.  When  the  local 
edge  vector,  as  defined  in  Table  1,  is  reversed. 
Equation  (6)  and  Equation  (7)  will  be  reversed  at  the 
same  time.  Therefore,  Equation  (6)  and  Equation  (7) 
should  be  multiplied  by  (-1)  if  the  local  edge  vector 
does  not  have  the  same  direction  as  the  global  edge 
direction. 

As  illustrated  in  Figure  4,  there  are  generally  four 
kinds  of  faces.  In  (a),  the  three  edges  of  the  face  are 
all  within  the  FEM  volume.  In  (b),  one  or  two  edges 
of  the  face  are  on  the  FEM/MoM  boundary.  In  (c), 
the  three  edges  of  the  face  are  all  on  the  FEM/MoM 
boundary  while  the  area  of  the  face  is  located  in  the 
FEM  volume.  In  (d),  the  three  edges  and  the  area  of 
the  face  are  all  on  the  FEM  and  MoM  boundary. 
Normally,  at  the  interface  between  higher-order  FEM 
elements  and  CN/LT  MoM  elements,  the  higher- 
order  terms,  e[  *  are  set  to  zero.  However,  for  the 

faces  of  type  (a),  (b)  and  (c),  the  complex  scalar 

rotating  along  the  edge  that  is  located  on  the 
FEM/MoM  boundary,  represents  fields  within  the 
FEM  volume  and  cannot  be  set  to  zero.  Allowing 
these  terms  to  have  a  non-zero  value  will  not  affect 
the  coupling  between  FEM  and  MoM,  since  their 
projection  on  the  boundaiy  is  equal  to  zero.  For  the 
faces  of  type  (d),  complex  scalars  £*  and 

Ef  corresponding  to  this  type  must  be  set  to  zero. 

IV.  NUMERICAL  RESULTS 

This  section  describes  three  examples  illustrating  the 
performance  of  the  hybrid  FEM/MoM  with  CT/LN 
and  LT/QN  FEM  basis  functions.  All  matrices  were 


solved  using  a  biconjugate  gradient  stabilized  solver 
[1 1].  A  750-MHz  Pentium  111  computer  was  used  to 
perform  the  computation. 


Figure  4.  Faces  in  a  tetrahedron. 

4.1  The  Scattered  Field  from  a  Sphere 

This  example  models  the  scattering  of  an 
electromagnetic  plane  wave  by  a  dielectric  sphere.  As 
shown  in  Figure  5,  the  radius  of  the  sphere  is  0.09  m 
and  the  relative  permittivity  is  4.5.  The  incident  wave 

propagates  in  the  +  z  direction.  The  wave  has 

amplitude  E0  and  is  polarized  in  the  x  direction, 

E(x,  y,  z,  /)  =  E0exeJM~fiz)  •  (23) 


WANG,  GUO,  HUBING:  HIGHER  ORDER  TANGETIAL  VECTOR  FEM/MOM  METHOD 


6 


Figure  5.  Scattered  field  from  a  dielectric  sphere. 


A  commercial  software  package  was  used  to 
discretize  the  FEM  volume  with  different  densities  to 
demonstrate  the  advantages  of  the  proposed  higher- 
order  TVFE’s.  The  MoM  boundary  was  chosen  to 
coincide  with  the  physical  boundary  of  the  dielectric 
sphere.  The  number  of  MoM  basis  functions  was 
fixed  during  the  whole  process.  For  validation,  results 
using  the  Mie  series  [15]  were  compared  to  the 
FEM/MoM  results.  In  Figure  6,  we  compare  results 
for  the  three-dimensional  bistatic  scattering  cross 
section  at  a  frequency  of  583  MHz.  The  Mie  series 
result  is  denoted  “Mie.”  For  a  mesh  with  a  small 
number  of  tetrahedra,  the  result  using  the  CT/LN 
TVFE  is  denoted  “CT/LN  TVFE  coarse,”  and  the 
result  using  the  LT/QN  TVFE  is  denoted  “LT/QN 
TVFE  coarse.”  For  a  mesh  with  a  larger  number  of 
tetrahedra,  the  result  using  the  CT/LN  TVFE  is 
denoted  “CT/LN  TVFE  dense.” 

In  Figure  6,  the  “CT/LN  TVFE  coarse”  result  is  seen 
to  compare  fairly  well  with  the  exact  Mie  series  result 
when  die  observation  angle  is  between  100  degrees 
and  180  degrees.  When  the  observation  angle  is 
below  100  degrees,  a  1-dB  discrepancy  can  be  seen 
because  die  mesh  is  relatively  coarse.  For  the  denser 
mesh,  die  “CT/LN  TVFE  dense”  result  shows  a 
significant  improvement.  By  keeping  the  original 
coarse  mesh  and  applying  the  LT/QN  basis  functions, 
the  “LT/QN  TVFE  coarse”  result  agrees  with  the 
exact  result  very  well.  Even  compared  with  the 
“CT/LN  TVFE  dense”  result,  the  “LT/QN  TVFE 
result”  is  closer  to  the  exact  result.  Table  2  presents 
relevant  parameters  for  die  three  results.  Improved 
accuracy  is  obtained  with  less  computer  resources 
using  LT/QN  FEM  basis  functions. 


Figure  6.  Bistatic  RCS  of  the  dielectric  sphere 
at  583  MHz. 


Table  2.  Comparison  between  the  results  in  Fig.  6 


FEM  Part 

FEM 

Unkno¬ 

wns 

MoM 

Unkno¬ 

wns 

Average 

Edge 

Length 

(mm) 

FEM 
Matrix 
Non  - 

zeros 

Solver 

Time 

(sec) 

CT/LN 

Coarse 

405 

346 

9.1 

2379 

28 

CT/LN 

Dense 

3266 

346 

4.3 

18336 

275 

LT/QN 

Coarse 

2430 

346 

9.1 

51052 

213 

4.2  Input  Impedance  of  a  Power  Bus  Structure 

This  example  models  a  printed  circuit  board  (PCB) 
power  bus  structure.  As  shown  in  Figure  7,  the  board 
dimensions  are  7.6  cm  x  5.1  cm  x  1.1  mm.  The  top 
and  bottom  planes  are  perfect  electric  conductors 
(PECs).  The  dielectric  between  the  PEC  layers  has  a 
relative  permittivity  of  3.81(l-j0.01).  The  MoM 
boundary  is  chosen  to  coincide  with  the  physical 
boundary  of  the  board.  A  source  is  identified  at  the 
location  shown  in  Figure  7. 


76  mm 


Source:  (28,  25.0)  and  (28,  25,  1.14) 


Figure  7.  A  PCB  power  bus  structure. 
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Unlike  the  previous  example,  the  fields  in  this 
configuration  are  relatively  uniform.  The  electric 
field  in  the  FEM  region  is  vertically  oriented  and 
constant  in  the  vertical  direction.  It  is  not  obvious 
that  a  higher-order  FEM  element  would  benefit  the 
analysis  of  this  configuration. 

The  FEM  uses  a  current  filament  on  tetrahedron 
edges  to  model  sources  located  within  the  FEM 
region  [16].  A  current  source  along  the  z-axis  can  be 
expressed  as, 


Jtol  =  lS(x  -  xf  )S(y  -yf)  z  (24) 

where  (xfi  yj)  specifies  its  position,  1  denotes  the 
electric  current  magnitude,  and  5(x)  is  the  Dirac  delta 
function.  The  contribution  to  vector  [g“"]  in  Equation 
(24)  is  simply, 

g“'  =  /  JJJ  :*{vi}S(x  -  xf  )S  (y  -  y  f )  dxdydz  •  (25) 

For  an  el  basis  function. 


For  an  e2  basis  function, 

For  fl  and  f2  basis  functions 
tu  =0 


(26) 


(27) 


(28) 


since  the  tangential  components  of  these  functions 
along  element  edges  are  zero. 


The  power  bus  structure  can  also  be  modeled 
analytically  as  a  cavity  with  two  PEC  and  four 
perfect  magnetic  conductor  (PMC)  walls.  The 
analytical  resonance  frequencies  are  given  as  follows 
(for  Pr=  10)  [17], 


where  a  and  b  are  the  length  and  width  of  the  cavity, 
respectively;  m  and  n  are  the  mode  indices;  c  is  the 
speed  of  light  in  free  space;  and  er  is  the  relative 
permittivity  of  the  material  in  the  cavity.  For  this 
power  bus  structure,  only  TMZ  modes  are  excited. 
From  Equation  (29),  the  TMZ  (1,0)  mode’s  resonance 
frequency  is  1011.1  MHz  and  foe  TMZ(2,0)  mode’s 
resonance  frequency  is  2022.3  MHz. 


(a)  TMZ(  1,0)  mode 


Figure  8.  Input  impedance  of  foe  power  bus  structure. 


In  Figure  8,  foe  computed  input  impedances  of  foe 
power  bus  structure  near  these  two  resonance 
frequencies  are  compared.  For  a  mesh  with  a  number 
of  tetrahedra  between  that  of  coarse  mesh  and  dense 
mesh,  foe  result  using  foe  CT/LN  TVFE  is  denoted 
“CT/LN  middle.”  The  “LT/QN  coarse”  and  “CT/LN 
coarse”  examples  employ  foe  same  mesh.  It  can  be 
seen  from  Figure  8  that  of  foe  four  cases,  the  “LT/QN 
coarse”  results  most  accurately  predict  the  resonance 
frequencies.  Table  3  presents  relevant  parameters  for 
foe  four  cases.  Once  again  the  LT/QN  coarse  mesh 
yields  more  accurate  results  with  fewer  computer 
resources  than  a  dense  CT/LN  mesh. 


Table  3.  Comparison  between  d 

ifferent  formulations 

Formulations 

Unknowns 

FEM 

non¬ 

zero 

Solver 

time 

(sec) 

CT/LN  Coarse 

467 

4097 

0.89 

CT/LN  Middle 

1162 

9940 

1.54 

CT/LN  Dense 

5041 

42039 

14.3 

LT/QN  Coarse 

4124 

87970 

6.8 
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In  [9],  the  inter-relationships  between  the  condition 
numbers  of  individual  elements  as  well  as  global 
matrices  based  on  various  interpolatory  and 
hierarchical  TVFE’s  were  studied  using  a  cavity 
resonator  example.  It  was  found  that  the  LT/QN 
TVFE  by  Andersen  and  Volakis  [7]  resulted  in  better 
conditioned  FEM  matrices  than  the  TVFE  by  Webb 
and  Forghani.  Since  resonant  structures  like  this 
power  bus  can  be  particularly  susceptible  to 
numerical  error,  this  structure  was  also  analyzed 
using  higher-order  elements  based  on  the  LT/QN 
elements  in  [7],  Figure  9  compares  the  condition 
numbers  of  the  global  FEM  matrices  generated  by 
this  example  based  on  Andersen’s  and  Webb’s 
LT/QN  basis  functions.  The  term  ‘norm’  denotes  the 
normalized  basis  functions  described  in  Appendix  A. 
Up  to  3  GHz,  the  LT/QN  TVFE  in  [7]  yields  slightly 
better  conditioned  matrices  than  the  TVFE  in  [5], 
Normalized  vector  basis  functions  yield  much 
smaller  condition  numbers  than  unnormalized  vector 
basis  functions. 


Figure  9.  Condition  numbers  of  the  global  FEM 
matrix  based  on  different  LT/QN  basis  functions. 


4.3  Input  Impedance  of  a  Microstrip  Structure 

Numerical  models  of  printed  circuit  board  (PCB) 
geometries  often  include  at  least  one  microstrip 
structure  (i.e.  a  trace  over  a  plane).  When  these 
structures  are  modeled  using  a  finite  element 
technique  with  CT/LN  basis  elements,  it  is  not 
uncommon  to  model  the  space  between  the  trace  and 
die  plane  with  a  single  layer  of  elements.  This 
approach  generally  yields  good  results  when  the  trace 
is  wide  (e.g.  the  power  bus  structure)  or  when  far- 
field  results  are  calculated.  However,  for  narrow 
traces  or  when  calculating  near-field  properties  (e.g. 
input  impedance  or  crosstalk),  a  single  layer  of 
elements  may  not  be  adequate  [16]. 


24.25  mm  !  5  mm  24  25  mm 


Figure  10.  The  geometry  of  a  microstrip  structure. 


Figure  10  shows  the  geometry  of  a  PCB  with  a  thin 
trace.  The  trace  width  is  only  slightly  greater  than  the 
trace  height.  The  board  is  made  of  a  dielectric  with 
£,=4.2.  The  trace  is  excited  by  a  current  source  at 

one  end,  and  is  terminated  by  a  47-ohm  resistor  at  the 
other  end.  The  MoM  boundary  is  chosen  to  coincide 
with  the  physical  boundary  of  the  board. 

The  FEM  code  models  load  impedances  zL  as 

dielectric  posts  on  tetrahedron  edges  [17],  Those 
posts  have  a  finite  conductivity  given  by 


where  /  is  its  length,  and  S  is  the  cross  sectional 
area.  If  the  load  is  treated  as  a  lumped  element,  its 
contribution  to  the  finite  element  matrix  is, 


[/,cl  ] = yr  HI  {w*‘  }•  {w"  r  s(*  -xLXy-  yL  )<k<fydz = \— 

zL 

(31) 

for  the  el  basis  function,  and, 

U'7 ]  =  yj.  J1/{W*J }•  {w'2 r -xLXy- yL )dxdydz 
=  7 -fa-Ltfdz 

(32) 

for  the  e2  basis  function,  and, 

k]=  0  (33) 

for  the  fl  and  f2  basis  functions.  The  electric  field 
lines  around  the  trace  are  illustrated  in  Figure  11 
[18].  Since  the  electric  field  around  the  trace  varies 
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dramatically,  the  coarse  mesh  used  to  divide  the  FEM 
volume  in  the  previous  power  bus  example  does  not 
work  here.  Figure  1 1  also  shows  two  meshes  for  a 
microstrip  geometry.  The  coarse  mesh  is  one  element 
tall  and  the  fine  mesh  is  two  elements  tall. 

Figure  12  shows  the  measured  and  calculated  results 
for  a  47-ohm  termination  up  to  1  GHz.  Note  that  the 
coarse  mesh  yields  a  poor  result  with  CT/LN 
elements  while  the  dense  mesh  results  are  close  to  the 
measured  results.  The  LT/QN  result  with  the  coarse 
mesh  also  yields  an  accurate  result. 

E  Field 
H  Field 


derivations  compared  to  other  hierarchical  higher- 
order  basis  functions.  The  properties  of  the  LT/QN 
basis  functions  are  discussed  and  compared  to 
traditional  CT/LN  basis  functions.  Three  examples 
demonstrate  that  higher-order  basis  functions  are 
capable  of  providing  more  accurate  results  with  a 
coarser  tetrahedral  mesh  and  less  computational 
resources.  The  condition  numbers  of  the  global  FEM 
matrices  derived  from  a  power  bus  structure  on  the 
basis  of  various  hierarchical  LT/QN  basis  functions 
are  compared.  It  is  confirmed  that  the  TVFE  by 
Andersen  results  in  somewhat  better  conditioned 
matrices  than  the  TVFE  by  Webb  and  Forghani. 
Also,  normalized  vector  basis  functions  are  observed 
to  result  in  much  smaller  condition  numbers  than 
unnormalized  vector  basis  functions. 


REFERENCES 

(a)  Electric  and  magnetic  field  lines 


0«a  sans 


(b)  Coarse  Mesh 
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APPENDIX  A 

FEM  analysis  requires  computation  of  two  matrices. 
These  two  matrices  are, 


£y  =  J  Vx  w,  •  Vx w;  dV  (Al) 

Fv  =  j  w(  •  w ;  dV  •  (A2) 

A.  CT/LN  TVFE 

(A3) 

B.  Webb's  unnormalized LT/QN  TVFE 
Edge  element 

w ?=l,lLdVLl2-Ll2VLll)  (A4) 

w :*« /,(£„?!„ +  £(1V£„)  /  =  l,-6  (A5) 

Face  element 

w,"  =  LnLl2VLl3  -  (A6) 

W,'!  =  LlXL,2VL,3-Ll2Ll3VL^  /  =  1,  — 4  (A7) 

B.  Webb ’s  normalized  LT/QN  TVFE 
Edge  element 

w =  L,]VLa-LaVL„  (A8) 

w^2  =  L„VL,2  +  L,2VLa  /  =  1,-6  (A9) 

Face  element 

=L,xLl2VLl3-LnLl3VLl2  (A  10) 

w /  =  1, — 4  (All) 

C.  Andersen's  unnormalized  LT/QN  TVFE 
Edge  element 

w  :'=ll(LnVL,2-Ll2VLll)  (A  12) 

<2  =  I. (4)  -  b,2 -  Z.,2  VL„ )  /  =  1,-6  (A  13) 
Face  element 

w  ?=LliLl2VLl3-LllLl3VLl2  (A14) 

'"l1-LALtlVLa-L,lLaVL,l  /  =  1, -4  (A  15) 

D.  Andersen  s  normalized  LT/QN  TVFE 
Edge  element 

W?  =i,,VLl2-i,jVi(l  (A  16) 

wj*  =(i„-i,!XW!-Wl)  /  =  l,"-6  (A  17) 
Face  element 

w  ?=LllLllVLl3-LllLl^L,2  (A  18) 

w r^L,2L,2Vll3-Ll2L^l,2  /  =  !,  -4  (A19) 
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Abstract.  The  calculation  of  integrals 
containing  the  free-space  Green’s  function  in 
electromagnetic  problems  is  difficult  to 
perform  with  great  accuracy.  Three 
approaches  to  the  calculation  are 
investigated.  The  inadequacy  of  the 
singularity-subtraction  method  is 
demonstrated.  The  Duffy  transform  is  shown 
to  provide  good  results  when  the  test-point  is 
on  the  surface  being  investigated.  A 
Maclaurin  series  expansion  with  integration 
prior  to  summation  is  shown  to  be  efficient 
and  reliable  both  on  and  off  the  surface  under 
study.  Solutions,  in  both  Cartesian  and 
cylindrical  coordinate  systems,  that  allow  the 
calculations  to  be  performed  to  a  pre-defined 
level  of  accuracy  are  presented. 

Introduction. 

The  magnetic  vector  potential,  or  MVP,  is  an 
important  quantity  that  appears  in  many 
electromagnetic  problems  that  involve 
evaluation  of  electric  and/or  magnetic  fields. 
For  example,  it  is  a  component  in  the 
definition  of  the  electric  field  integral  equation, 
EFIE  [1,  pi  7].  In  this  context  it  may  exist  in 
its  basic  form  or  it  may  be  subject  to 
differentiation.  Its  use  is  at  its  most  basic 
when  used  in  the  solution  of  Hallen’s  integral 
equation  for  a  cylindrical  dipole.  When  used 
for  solving  Pocklington’s  integro-differential 
equation  for  the  same  dipole  the  second 
differential  of  the  MVP  must  be  considered. 
Derivatives  of  the  MVP  are  also  derivatives  of 
the  Green’s  function  contained  within  the 
definition  of  the  MVP.  Because  the  three- 
dimensional  Green’s  function  contains  a 
singularity,  it  is  preferable  to  keep  the  order  of 


differentiation  to  a  minimum,  preferably  zero. 
When  it  cannot  be  kept  to  zero,  then  one  of 
two  actions  are  generally  undertaken.  Either 
the  derivatives  must  be  transferred  to  the 
basis/testing  function  used  in  the  solution  of 
the  particular  EFIE  under  investigation  or  one 
of  the  special  formulations  that  have  been 
developed  to  accommodate  the  differentiation 
[2]  [3]  must  be  considered.  Even  when  one 
examines  the  evaluation  of  just  the  non- 
differentiated  form  of  the  MVP  one  has  certain 
numerical  difficulties  to  face.  These 
difficulties  are  addressed  in  this  report. 
Accurate  evaluation  of  the  MVP  is  gaining 
importance  as  the  use  of  higher  and  higher 
order  basis  functions  is  considered.  Also,  as 
we  shall  see,  evaluation  of  the  MVP  when  the 
test-point  is  located  a  short  distance  from  the 
test  surface  is  a  requirement  that,  while  of 
interest  in  many  applications,  is  handled 
poorly  by  current  techniques. 

As  the  title  of  this  paper  suggests,  the  over¬ 
riding  issue  here  is  one  of  solution  accuracy. 
When  calculating  entries  in  a  matrix,  Z ,  and 
then  solving  the  corresponding  matrix 
equation,  Z7  =  V ,  Miller  [4]  has  shown  that 
the  solution  error,  |rf/|/|/|,  is  comparable  to 

the  product  of  the  error  in  the  terms  in  the 
matrix  and  the  condition  number  of  the  matrix. 
For  example,  if  the  error  of  the  matrix  terms, 

||</Z||/||Z||,  is  10-6  and  the  condition  number  is 

104 ,  the  error  of  the  solution  may  be  no  better 

than  10"2  -  if  nothing  else  introduces  further 
errors.  The  resulting  accuracy  is  at  the  lower 
limit  of  usefulness.  It  can  be  improved  by 
either  reducing  the  condition  number  of  the 
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matrix  or  by  reducing  the  error  in  the 
calculation  of  the  matrix  entries.  This  work  is 
about  the  latter. 

After  setting  out  some  basic  definitions,  this 
report  will  investigate  the  difficulties 
associated  with  the  most  widely  used 
approach  to  the  evaluation  of  the  MVP.  The 
second  section  will  examine  the  use  of  the 
Duffy  [5]  transform  which  was  originally 
conceived  to  address  issues  arising  in  the 
evaluation  of  integrals  such  as  found  in  the 
MVP.  The  third  section  will  show  how  the 
MVP  can  be  dealt  with  in  a  manner  that  is  not 
only  rigorous  but  also  efficient.  This  will  be 
followed  by  a  discussion  of  testing, 
concluding  with  a  statement  of  key  findings. 

Definitions. 


The  three-dimensional  MVP  is  defined  as: 


When  examining  currents  on  surfaces,  this 
definition  reduces  to: 


(2) 


In  the  above,  r  denotes  the  position  vector  to 
the  test/observation  point,  and  r'  is  the 
position  vector  to  the  surface  under  study. 
The  scalar  component(s),  /4. ,  of  the  surface 


current  will  be  represented  by  polynomials,  in 
u,  in  the  form:  Is  =  . 

i=0 


Machine  precision  will  be  referred  to  often  in 
this  study.  By  machine  precision  we  will  be 
referring  to  machine  epsilon,  e ,  which  is  the 
gap  between  1  and  the  next  larger  floating 
point  number  [6,  pi 4],  e  =  2~(p~])  where  p 
is  the  precision  of  the  machine  in  bits. 
Machine  precision  is: 


log10(*)  =  ~(P  ~  l)log10(2)  (3) 

Results  for  Compaq  Fortran  on  an  Alpha 
processor  are  shown  in  Table  I. 


Relative  error,  Rel.  error,  is  defined  as: 


Rel.  error  =  '/a/T  M- 


fref  will  be  defined  each  time  that  relative 

error  is  discussed.  When  evaluating  integrals 
one  frequently  compares  the  result,  /„ , 


obtained  in  the  most  recent  evaluation,  to  the 
result,  obtained  in  the  prior  evaluation. 

This  is  more  accurately  defined  as 
convergence  rate  and: 


Convergence  Rate  = 


(5) 


Singularity  Subtraction. 


The  free-space  Green’s  function  is  defined  as 
G(R )  =  e~JkR/R  where  R  is  the  distance 
between  the  source  and  the  observation  or 
test  point.  The  mathematical  definition  of  R 
is  specific  to  the  coordinate  system  in  use  and 
will  be  elaborated  on  later.  For  this  section  of 

the  report  we  define  R  =  V«2  +  S2  where  u 
is  an  independent  variable  and  S  will  assume 
various  fixed  values.  G(R)  is  split  into  two 

parts  G(R)  =  Gr(R )  +  GS(R),  where 
Gr(i?)is  the  non-singular  part  and  GS(R)  is 
the  singular  part  [7],  Specifically,  we  have: 


Gr(R )  =  e~JkR/R  -  ]/R0 

(6a) 

GS(R)  =  \/R0 

(6b) 

GS(R)  is  developed  from  a  Taylor  series 
expansion  of  G(R)  and  when  possible  it  is 
evaluated  analytically  [8].  An  example  of  a 
solution  in  the  Cartesian  coordinate  system  is 
given  in  [1,  p420].  In  the  event  an  analytical 
solution  is  not  available,  a  recent  discussion 
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of  numerical  methods  is  available  in  [9].  The 
focus  in  the  following  is  on  Gr(R ) . 


For  purposes  of  this  immediate  analysis  we 
express  Gr(R)  in  terms  of  its  real  and 
imaginary  components: 

Gr(R)  =  (cos (kR)  -  1  )/R  -  jsm(kR)/R  (7) 

These  two  components,  when  R  =  |«|  (i.e. 

8  =  0)  are  plotted  in  Figure  1.  Both 
components  are  finite  throughout  the  range. 
However,  the  real  component  is  obviously  not 
‘smooth’  at  u  =  0.0 . 

The  results  for  integrating  the  real  component 
of  Gr(R)  with  Gauss-Legendre  quadrature 
are  plotted  in  Figure  2a,  where  d  is  the  same 
as  8  in  the  text.  The  reference  values  were 
calculated  using  the  series  expansion  method 
described  later.  It  is  observed  that  the 

integration  convergence  becomes  worse  as 
S  -»  0,  a  finding  which  is  somewhat 
unexpected,  counter-intuitive  and 

disconcerting.  A  similar  observation  is 

implied  in  Figure  2  in  [10]. 

When  evaluating  integrals  numerically, 
particularly  close  to  the  source,  it  is  important 
to  remember  that  integration  rules  generally 
exhibit  an  error  that  is  proportional  to  the 
derivatives  of  the  function  being  integrated. 
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In  the  case  of:  /  =  (cos(kR)  -  1  )/R  (8a) 


From  this  last  equation  we  see  that  as: 

which  provides  an  explanation  for  the 
behavior  of  the  curves  in  Figure  2a.  To 
illustrate  that  it  is,  in  fact,  the  presence  of  the 
discontinuity  at  u  =  0  that  creates  the 
problems,  the  lower  limit  of  integration  was 
moved  from  u  =  Otow  =  0.01.  The  results 
of  this  change  are  shown  in  Figure  2b  and 
clearly  demonstrate  that  the  problem  has 
been  significantly  mitigated. 

For  the  special  case  of  8  =  0 ,  we  find  that 
although  there  is  a  jump  in  the  values  of 

df  d 2  f 

— ,  — y  and  higher  derivatives  when 

du  du 

u  moves  from  0_  to  0+  they  are  finite  in  that 
region  nevertheless,  and  so  the  integration 
rules  hold  for  this  special  case. 

For  double/triple  integrals  8  is  introduced  into 
the  inner  integral  by  the  outer  integral  defining 
the  MVP  and  thus  8  may  become  arbitrarily 


*i »  e~jkr  &c  e~jkr  ylK^-jkr  - 

jj - dxdy  =  ]dx  J - dy  +  J dy  J - dx  where  r  =  yjx2  +  y2  +  z2  and  AT  =  yj  jc, 

00  r  oor  0  0  r 

Substituting  y  =  uKx  and  x  =  vy/K  in  the  inner  integrals,  we  arrive  at: 

-  Yfrju _ Ct _ _  +  UL — £*l - - 

0  0  V«2  +  {MKf  +  {zJKxf  1  Uv2  +K2  +  (z0K  /  y)2 


nil  0/K  -fir 


|  jzp— — dddz  =  j  dO  J  zp- — dz  +  jzpdz  j——d0 


where  r  =  Jz2  +  A p\  +  4 a(a  +  Ap0)sin2^,  K= — 

2z, 

Substituting  z  =  uO/K  and  ^  =  vKz  in  the  inner  integrals  we  obtain: 

i*  -  1m  ,  - ,  + 

0  oJu2  +  4  a(a  +  &p0)(K  sin  d/0)2  +  (KApJd)2 


Jzpdzj- 


oJ(l/Ky  +  4  a(a  +  Ap0)(sm(vKz)  I  Kzf  +  (A pJKz)1 


Equation  Set  1 .  The  MVP  equation,  expressed  in  two  different  coordinate  systems, 
transformed,  by  the  Duffy  method,  to  remove  the  singularity  at  the  origin. 
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small,  but  non-zero,  leading  to  these 
problems.  Consequently,  when  evaluating 
the  real  part  of  Gr(R ) ,  it  would  appear  that  all 
quadrature  rules,  when  applied  directly,  are 
doomed  to  fail  as  even  the  simple  trapezoid 
rule  requires  that  the  second  derivative  of  the 
integrand  be  well  behaved.  With  this 
conclusion  it  is  advisable  to  seek  alternative 
methods. 

Singularity  Removal  by  Transformation. 

In  1982,  Duffy  [5]  proposed  a  method  that, 
through  a  change  of  variable,  causes  the 
removal  of  the  singularity  in  the  integrand  of 
two  and  three  dimensional  integrals.  His 
method  is  presented  here  first  with  a  constant 
current,  for  simplicity,  in  the  Cartesian 
coordinate  system.  Currents  of  polynomial 
form  up ,  are  included  with  the  discussion  of 
the  cylindrical  coordinate  system. 

The  formalism  for  each  of  the  two  coordinate 
systems  is  shown  in  Equation  Set  1 .  It  allows 
for  the  test/observer  point  to  be  offset  from 
the  surface  -  by  amount  z0  in  the  Cartesian 

system  and  A p0  in  the  cylindrical  system. 

When  the  offset  is  zero,  the  formulae  in  the 
Cartesian  system  clearly  show  that  the 
singularity  has  been  removed.  Furthermore, 
in  this  case,  the  denominator  in  the  integrand 
is  not  dependant  on  the  variable  associated 
with  the  outer  loop  other  than  through  the 
value  of  K,  which  is  fixed.  Consequently, 
the  derivatives  of  the  integrand  are  all  well 
behaved  and  one  can  expect  that  quadrature 
integration  will  work  well,  and  indeed  it  does. 


Precision 

€ 

Outer 

Integral 

Inner 

Integral 

Single 

1.2E-07 

7 

7 

Double 

2.3D-16 

7 

13 

Quad 

2.0D-34 

13 

37 

Table  1.  T 
one  of  the 
domain  usin 

he  number  of  terms  required  in 
double  integrals  for  a  square 
g  the  Duffy  transform. 

For  example,  using  Gauss-Legendre  for  both 
the  inner  and  outer  integrals,  the  numbers  of 
terms  necessary  for  computing  the  double 
integral,  to  a  precision  of  2s ,  on  a  flat  surface 
are  shown  in  Table  1 .  The  dimensions  of  the 
cell  were  0.0  <  x  <  0.1,  0.0  <  y  <  0.1. 

Because  of  symmetry  in  this  example,  the 
numbers  of  terms,  required  in  the  two  double 
integrals,  are  the  same.  All  the  integrals 
terminate  when  the  convergence  rate  falls 
below  the  precision  level,  2s . 

When  examining  a  cylindrical  case  we  look  at 
a  cell  width  of  0.1  wavelengths  on  a  cylinder 
of  radius  0.007  wavelengths.  In  addition, 
polynomial  representations  of  the  current,  to 
the  degree  p,  are  incorporated.  Again,  the 
examination  takes  place  on  the  surface,  so 
that  Ap0  =  0.  For  this  case  we  track  the 

number  of  terms  in  the  outer(out)  and 
inner(in)  integrals  for  the  two  transformed 
double  integrals.  We  identify  these  as  u-out, 
u-in,  v-out  and  v-in.  The  results,  calculated  in 
double  precision  to  a  precision  level  of  2s, 
are  shown  in  Table  II.  They  indicate  no 
dependence  between  the  degree,  p,  and  the 
number  of  terms. 


JL_ 

u-out 

u-in 

v-out 

v-in 

Total 

0 

11 

36 

8 

10 

476 

1 

9 

38 

9 

10 

432 

2 

10 

36 

10 

9 

450 

3 

13 

34 

9 

10 

532 

4 

12 

29 

10 

9 

438 

5 

11 

29 

10 

9 

409 

Table  II.  1 

fhe  number  of  terms  required  in  1 

each  of  the  integrals  for  the  cylindrical 
surface  using  the  Duffy  transform. 

The  good  performance  of  the  Duffy  transform 
when  the  offset  is  zero  does  not  follow 
through  when  the  offset  is  some  finite  value. 
This  is  revealed  when  the  data  of  Table  III  is 
examined.  This  table  shows  the  total  number 
of  iterations  needed  by  the  integrals  as  a 
function  of  the  value  of  the  offset,  A p0 .  The 

dimensions  on  a  cylindrical  surface  are  the 
same  as  used  in  Table  II.  A  value  of  p  =  0 
was  used. 
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Offset 

Duffy 

G-L 

0.0 

476 

N/A 

0.0001 

12280 

4356 

0.001 

4504 

1860 

0.01 

1370 

952 

0.1 

311 

143 

Table  III.  The  effect  of  the  offset  value  on 
the  number  of  integration  terms  for  the 
cylindrical  surface  using  the  Duffy  transform 
and  Gauss-Legendre  applied  directly  to  the 
double  integral. _ 

For  comparative  purposes  Table  III  also 
shows  the  results  when  the  MVP  is  integrated 
directly.  This  is  possible  when  the  offset  is 
finite.  We  conclude  that  for  test/observation 
points  off  the  surface,  the  Duffy  transform  is 
unacceptably  inefficient.  Nevertheless,  when 
the  test/observation  point  is  on  the  surface 
the  transform  offers  a  method  that  provides 
rigorous  convergence  in  the  integrals 
associated  with  the  MVP. 

A  Series  Expansion  for  G(R ) . 

The  Maclaurin  expansion  for  the  real  and 
imaginary  components  of  the  Green’s  function 
are: 

cos (kR)/R  =  1 /R  -  k2R/ 2!  +  ^4i?3/4!..(1 1  a) 

sm(kR)/R  =  k  -  +  k5RA/5\...  (11b) 

The  method  proposed  integrates  the 
expansions  in  (11a)  and  (11b),  term  by  term, 
until  the  ratio  of  the  last  term  evaluated  to  the 
largest  term  evaluated  is  less  than  machine- 
precision.  In  this  way  it  is  possible  to 
develop  analytical  terms  for  the  inner  integral. 
The  terms  for  the  expansion  in  (11a)  are 
shown  in  the  series  (12a)-(12d)  and  the 
terms  for  the  expansion  of  (11b)  are  shown  in 
(13a)-(13d).  Each  of  these  integrals  is  exact 
for  a  given  value  of  8 .  The  formulae  are 
shown  with  a  lower  limit  of  u  =  0.0 ,  but  this 
is  done  for  convenience  only. 


The  Green’s  function  is  rarely  evaluated  on  a 
stand-alone  basis;  rather  it  is  evaluated  in 
conjunctions  with  a  representation  of  the 
current  on  the  surface  being  studied.  We  will 
adopt  the  polynomial  summation  defined 
earlier.  The  integrals  of  interest  are: 

^  -  y^du  (16a) 

Gta  =  j (i6b) 

i  R 

The  case  of  p  =  0  has  already  been 
presented  in  equations  (12)  and  (13).  The 
case  for  p  =  1  is  shown  in  equation  (14)  and 
for  p  >  2  the  relevant  equation  is  (15).  This 
last  equation  is  applicable  to  the 
computations  of  both  GRe  and  Gta .  Thus, 
once  the  terms  for  p  =  0  and  p  =  1  have 
been  evaluated,  the  evaluations  for  higher 
values  of  p  are  straightforward. 

The  Series  Expansion,  Integration  and 
Summation,  SEIS,  process  described  above 
was  tested  in  several  ways.  The  first  test 
revisits  the  calculations  performed  for  use  in 
Figure  2a.  We  examine  the  effect  of  varying 
8  on  the  number  of  terms  required  to  achieve 
convergence  in  the  summations  of  the  series. 
The  results  of  such  calculations  on  GRe  for 
p  =  0 ,  shown  in  Figure  3,  demonstrate  that 

the  value  of  8  has  little  impact  in  this  context 
and  thus  we  conclude  that  a  major  goal  of  the 
present  work  has  been  achieved. 

The  second  series  of  tests  performed 
involved  the  inclusion  of  basis  functions  as 
discussed  earlier.  Shown  below,  in  Equation 
Set  III,  are  analytical  expressions  for  some 
inner  integrals.  Using  these  expressions,  the 
accuracy  of  the  present  method  can  be 
examined  for  the  inner  integral  when  p  is  odd. 
As  a  practical  note,  the  calculation  of  these 
inner  integrals  (17)  to  high  accuracy  required 
the  use  of  extremely  high  precision  software. 
That  used  here  was  developed  by  Bailey  [1 1], 
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=  log (b  +  sfb2  +  8s)  -  log(|4 

o  ° 

b  1  r  h  ri 

^Rdu  =  —  bRb  +  S2  J—  where  Rh  - 

o  ^  _  0*. 

b  ,  r  h 

J R3du  =  -  bRl  +  3 S2  \Rdu 


=  sib 2  + 


/>  -j  b 

f R2m-'du  =  —  bRlm-'  +  (2w  -  1)£2  {R2m-3du  where  m  >  1 

0  M  0 


(12a) 


(12b) 


(12c) 


(1 2d) 


o 

jdu  =  b 
o 

A  A  ,3 

J  R2du  =  J(w2  +  S2)du  =  —  +  <?2& 

0  0  3 

JtfVw  =  J(M2  +  S2)2du  =  —  +  2  82—  +  b 


(13a) 


(13b) 


(13c) 


fl>2n  ,  f,  2  jOv,  ,  6""' 

=  (w  +  S')  du  = - + - . 

i  j  2n  +  1  2(2 n  - 1) 


+  82nb 


len  p  =  1 


(13d) 


»2  nm+l  ““"2 

[w/T"1**/  = -  where  0  <  m 

•  m  +  1 

Wi  «=«, 

len  p  >  2 ,  a  recurrence  formula  can  be  derived  which  takes  the  form: 

\upRmdu  =  - - - -  up-'Rm+2  “=“a  -(p-  1  )82\up-2Rmdu  where 

J  p  +  m  +  1[  "=«.  J 


Rmdu  where  -  1  <  m  (15) 


Equation  Set  II.  The  basic  equations  for  the  term-by-term  integration  of  a  Maclaurin 
series  expansion  of  the  Green’s  function  and  its  product  with  a  polynomial. 


The  results  for  GRe  appear  in  Figure  4. 
Results  for  Glm  and  GReare  not  visually 
distinguishable  and  hence  only  the  results 
for  GRe  are  shown.  The  plots  clearly 

show  that,  at  most,  5,  8,  or  14  terms,  for 
single,  double  or  quad  precision 
respectively,  are  needed  in  the  series 
expansion  for  this  integration  range.  The 
plots  also  clearly  show  that  the  relative 


error  is  not  dependant  on  p,  the  exponent 
in  the  basis  function  used  to  represent  the 
current.  The  results  presented  in  Figure 
5,  for  the  even  values  of  p,  are  referenced 
with  respect  to  their  own  machine 
precision  limited  values.  Their  behavior 
is  similar  to  the  results  for  odd  values  of  p, 
which  are  referenced  to  analytical  values. 
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%upP~jkR  i 

onsider  I  =  f - du  where  R  =  Vw2  +  82 .  With  the  substitution  dR  = 

p  J  J? 


we  get  lp  =  J  up~'e'JkRdR.  When  p  is  odd,  this  integral  can  be  solved  analytically. 
For  example: 

i-  i  j  e  'M  P 

p  =  1’  4  =  ~r~ 


(17a) 


P  =  3;  / 


=  e~M  u2\  —  \  -  2R\  —  ]  +  2| — 


P  =  5;  /.  = 


w4  -  4  +  (£2  -  37?  )  -  +  67?l  f-  -  6 


(17b) 


(17c) 


Equation  Set  III.  Examples  of  analytical  solutions  for  the  integral  for  odd  p. 

The  surface  is  that  of  a  cylinder  of  radius  a, 
Application  to  the  Magnetic  Vector  and  with  A p0  >  -a  the  test  point  is  at 


Potential. 


(a  +  A/70,0,z0). 


In  this  instance 


The  application  of  the  above  procedure  to 
double  integrals  is  straightforward.  The  inner 
integral  is  computed  as  above  and  then  a 
quadrature  integration  formula  is  applied  to 
the  outer  integral. 


Cartesian  coordinates.  A  double  integral  of 
interest  is  given  by: 


where 


R  =  J(x  -  x0)2  +(y  -  _y0)2  +  4 


The  surface  is  in  the  x-y  plane  and  the  test 
point  is  at  (x0,>>0,z0) .  In  terms  of  the  inner 

integral  we  replace  (x  -  x0)  by  u , 


8 2  =  (A/?0)2  +  4 a{a  +  Ap0)s\n2  <p  and 
u  =  (z  -  z0)  and  again  the  limits  of 
integration  are  appropriately  adjusted. 

The  results  for  the  calculation  of  the  MVP  for 
a  section  of  a  cylindrical  dipole  with  values  of 
a  =  0.007>l  and  0.0  <  z  <  0.1  are  shown 
in  Figures  8a  and  8b.  Two  quadrature 
methods  were  investigated  -  the  Gauss- 
Legendre  method  and  the  Linlog  method  [12]. 
The  reason  for  choosing  the  latter  method  is 
that  the  series  containing  even  p  always 
contains  log  terms  in  its  real  component. 
Linlog  was  designed  specifically  to  integrate 
functions  that  contain  polynomials  and 
logarithmic  terms. 


(O'  -  ^o)2  +  zl)  by  8 2  and  adjust  the 
integration  limits  appropriately. 

Cylindrical  coordinates.  Here,  the  double 
integral  is  given  by: 

1  2r  "r  e~jkR 

4  v  =  —  J^J(z  -  Z«Y  where  R  = 
V(z  -  ~zo)2  +  (APo )2  +  4«(«  +  A/?0)sin2^ 


Both  sets  of  calculations  were  performed  in 
quad  precision.  The  reference  values  were 
calculated  using  42  terms  with  the  respective 
quadrature  methods.  The  superiority  of  the 
Linlog  approach,  when  applied  for  even 
powers  of  p,  is  clearly  visible.  The  relative 
error  is  seen  to  reach  a  level  of  approximately 
-21 .4  and  then  remains  constant.  The  nodes 
and  weights,  as  originally  reported,  are  only 
known  to  20  digits,  hence  the  observation  is 
hardly  surprising.  This,  then,  is  the  bound  on 
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the  accuracy  to  which  this  particular  MVP  can 
be  calculated  using  today’s  tools. 

It  remains  to  examine  the  effect  of  the  offset 
value  as  was  done  for  the  Duffy  transform 
and  reported  in  Table  III.  Again  the  work  is 
done  in  double  precision  and  is  for  p  -  0. 
The  ensuing  results  are  shown  in  Table  IV. 
Compared  with  the  results  in  Table  III,  it  is 
clear  that  the  value  of  the  offset  has  little 
effect  on  the  number  of  terms  needed  to 
achieve  a  relative  accuracy  equal  to  the 
machine  precision. 


Offset 

Outer 

Series 

Total 

0.0 

14 

8 

112 

0.1 

17 

9 

153 

0.01 

23 

7 

161 

0.001 

24 

8 

192 

0.0001 

24 

7 

168 

Table  IV. 

The  effect  of  the  offset  values  on  1 

the  number  of  integration  terms  for  the 

cylindrical  surface  using  the  term-by-term 
integration  of  a  Maclaurin  series. 

Comparison  Between  Duffy  and  SEIS. 


The  efficiency  of  the  calculation  of  the  MVP 
by  the  two  methods  -  the  Duffy  transform  and 
the  SEIS  method  -  was  investigated  for  the 
cylindrical  case  already  discussed.  In  the 
case  of  the  Duffy  transform  the  number  of 
function  evaluations  of  both  inner  integrals 
was  counted.  In  the  case  of  the  SEIS  the 
count  was  the  product  of  the  number  of  nodes 
in  the  outer  integral  and  the  number  of  terms 
in  the  series  expansion.  The  results  are 
shown  in  Figure  9.  The  reference  line  is 
located  at  jra.  It  appears  that  both  methods 
are  most  efficient  when  the  aspect  ratio  of  the 
cell  under  consideration  is  approximately  1:1. 
The  Duffy  transform  is  particularly  susceptible 
to  this  phenomenon.  In  the  case  of  the  SEIS 
approach,  at  the  high  end  of  the  z  range,  the 
number  increases  as  the  value  of  z  increases 
-  due  to  the  need  for  more  terms  in  the  series 
expansion.  At  the  low  end  of  the  range,  the 
number  of  terms  needed  in  the  series 
expansion  falls  off  -  but  the  number  of  nodes 
needed  in  the  Linlog  integration  increases. 


z  dimension,  wavelengths,  on  the  cylinder 


Figure  9.  The  effect  of  the  z  dimension  on 
the  number  of  'evaluations'  needed  to 
compute  the  MVP  using  the  Duffy  transform 

and  the  SEIS  methods. 

n  this  example,  the  number  of  integration 
nodes  needed  for  z=0.001  was  39. 

Key  Findings 

1)  It  was  shown  that  conventional 
numerical  methods  give  misleading  results 
when  integrating  the  Green’s  function. 
Consider  the  results  of  Figure  2a.  When 
8  =  0.0001  the  relative  error  changes  very 
little  as  the  number  of  integration  terms  is 
increased  until  very  large  numbers  of  terms 
are  employed.  The  slow  improvement  in  the 
error  curve  would  normally  be  interpreted  as 
convergence  -  leading  to  an  inaccurate 
evaluation.  It  was  shown  that  this  behavior  is 
a  direct  consequence  of  the  derivatives  in  the 
neighborhood  of  u  =  0 . 

2)  The  results  presented  in  Figures  2a 
and  6  emphatically  illustrate  the  poor 
performance  of  Gauss-Legendre  methods  for 
evaluating  the  real  component  of  any  of  the 
integrals  studied.  This  finding  is  applicable  to 
all  quadrature  methods  that  are  applied 
directly  to  this  class  of  problem. 
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3)  The  Duffy  transform  provides  a 
reliable  method  for  computing  the  MVP  when 
the  test  point  is  on  the  surface.  As 
implemented  here,  it  is  not  suitable  for  use 
when  the  test  point  is  off  the  surface.  It  has 
the  advantage  that  the  integrals  can  be 
evaluated  using  standard  integration 
techniques  such  as  Gauss-Legendre. 

4)  The  use  of  a  Maclaurin  expansion  of 
the  Green’s  function,  followed  by  term  by  term 
integration  and  careful  summation  provides  a 
stable  means  for  calculating  both  the  real  and 
imaginary  components  of  the  function.  The 
method  is  efficient  and  can  be  used  both  on 
and  off  the  surface  being  examined. 

5)  Analytical  solutions  for  the  integral  of 
the  Green’s  function  and  its  product  with 
polynomial  representations,  of  odd  degree,  of 
the  surface  current  have  been  presented. 
These  solutions  provide  a  method  for 
evaluating  both  the  convergence  and  the 
accuracy  of  the  series-expansion-integration- 
summation  approach. 

6)  The  analytical  results  presented  for 
current  representations  of  odd  degree  would 
appear  to  offer  an  accurate  and  efficient 
approach  to  the  evaluation  of  those  integrals. 
However,  it  was  found  that  rounding  errors 
seriously  degraded  the  accuracy  of  such 
calculations  when  the  range  of  integration 
was  small  and  such  an  approach  should  be 
avoided  unless  high  precision  software  is 
employed. 

7)  The  algorithm  used  for  the  outer 
integral,  when  using  the  series  expansion- 
integration-summation  approach  must 
recognize  the  presence  of  the  logarithmic 
terms  in  the  series  expansion  when  p  is 
even.  This  means  using  the  Linlog  method 
[12]  for  this  particular  integration. 

Final  Remarks 

The  calculation  of  integrals  associated  with 
the  magnetic  vector  potential  has  been 
examined  in  depth.  The  integration  of  the 
Green’s  function  should  not  be  attempted  with 


quadrature  methods,  unless  some  suitable 
transformation  to  remove  the  singularity  has 
been  undertaken.  An  example  of  the  latter  is 
the  transform  due  to  Duffy,  and  this  is  quite 
suitable  when  the  test  point  is  on  the  surface. 
For  all-round  performance,  it  is  proposed  that 
the  inner  integral,  that  includes  the  Green’s 
function,  be  evaluated  by  means  of  the 
integration  of  each  term  of  a  Maclaurin 
expansion.  The  outer  integral  can  then  be 
evaluated  using  the  Linlog  rule.  In  all  cases, 
the  integration  can  be  taken  to  the  precision 
of  the  machine/compiler  (single,  double  or 
extended/quad),  except  that  the  Linlog 
nodes/weights  currently  limit  the  relative  error 
to  approximately  20  digits.  The  integration  of 
the  Maclaurin  series  prior  to  summation 
provides  a  method  that  is  efficient  and 
accurate. 

References. 

[1]  A.  F.  Peterson,  S.  L.  Ray,  R.  Mittra, 
Computational  Methods  for  Electromagnetics, 
IEEE  Press,  1998. 

[2]  D.  B.  Miron,  “The  Singular  Integral 
Problem  in  Surfaces”,  IEEE  Trans.  Ant.  & 
Prop.,  Vol.  AP-31,  No.  3,  pp.  507-509,  May 
1983. 

[3]  A.  D.  Yaghijan,  “Electric  Dyadic 
Green’s  Functions  in  the  Source  Region”, 
Proc.  IEEE,  Vol.  68,  No.  2,  pp.  248-263,  Feb. 
1980. 

[4]  E.  K.  Miller,  “A  Computational  Study  of 
the  Effect  of  Matrix  Size  and  Type,  Condition 
Number,  Coefficient  Accuracy  and 
Computation  Precision  on  Matrix-Solution 
Accuracy”,  Inti.  Symp.  On  Ant.  &  Prop.  Digest, 
Newport  Beach,  CA,  pp.  1020-1023,  June 
1995. 

[5]  M.  G.  Duffy,  “Quadrature  Over  a 
Pyramid  or  Cube  of  Integrands  with  a 
Singularity  at  a  Vertex”,  SIAM  J.  Numer. 
Anal.,  Vol.  19,  No.  6,  pp.  1260-1262,  Dec. 
1982. 

[6]  M.  L.  Overton,  Numerical  Computing 
with  IEEE  Floating  Point  Arithmetic,  SIAM, 
2001. 

[7]  D.  R.  Wilton,  C.  M.  Butler,  “Effective 
Methods  for  Solving  Integral  and  Integra- 


BIBBY,  PETERSON:  CALCULATION  OF  MAGNETIC  VECTOR  POTENTIAL 


22 


Differential  Equations",  Electromagnetics,  Vol. 
1,pp.  289-308,  1981. 

[8]  J.  M.  Song,  W.  C.  Chew,  “Moment 
Method  Solutions  Using  Parametric 
Geometry",  J.  of  Electromagnetic  Waves  and 
Appl.,  Vol.  9,  pp.  71-83,  1995. 

[9]  A.  Herschlein,  J.  v.  Hagen,  W. 
Wiesbeck,  “Methods  for  the  Evaluation  of 
Regular,  Weakly  Singular  and  Strongly 
Singular  Surface  Reaction  Integrals  Arising  in 
Method  of  Moments”,  ACES  Journal,  Vol.  15, 
No.  1,  pp.  63-73,  March  2002. 

[10]  M.  Gimersky,  S.  Amari,  J.  Bornemann, 
“Numerical  Evaluation  of  the  Two- 
Dimensional  Generalized  Exponential 
Integral”,  IEEE  Trans.  Ant.  &  Prop.,  Vol.  AP- 
44,  No.  10,  pp.  1422-1425,  1996. 

[11]  D.  H.  Bailey,  “A  Portable  High 
Performance  Multi-precision  Package",  RNR 
Technical  Report  RNR-90-022,  May  18,  1993. 
See  also  www.nersc.gov/~dhbailey/mpdist/. 

[12]  J.  H.  Ma,  V.  Rohklin,  S.  Wandzura, 
“Generalized  Gaussian  Quadrature  Rules  for 
Systems  of  Arbitrary  Functions",  SIAM  J. 
Numer.  Anal.,  Vol.  33,  pp.  971-996,  June 
1996. 

Addendum:  Procedural  Considerations. 

In  order  to  assure  the  best  possible  accuracy 
when  evaluating  series  such  as  those  in 
equations  (12)  -  (15),  several  “good  practice” 
issues  need  to  be  followed. 

1)  Terms  that  are  to  be  added  need  to  be 
stored  separately  from  those  that  are  to  be 
subtracted.  Thus  the  values  associated  with 
the  upper  limit  in  an  integration  formula  must 
be  separated  from  those  associated  with  the 
lower  limit  of  integration.  A  similar  separation 
should  be  maintained  when  implementing  the 
Maclaurin  series,  noting  that  this  involves  the 
additional  complication  of  a  series  with 
alternating  signs. 

2)  When  evaluating  the  terms  associated  with 
the  real  part  of  the  overall  integral,  negative 
values  may  occur  due  to  the  presence  of  the 
log  term.  These  should  be  identified  and 
stored  appropriately. 

3)  Terms  should  be  added  by  starting  with  the 
smallest  and  proceeding  to  the  largest.  To 


this  end,  the  two  sets  of  terms  need  to  be 
sorted  in  ascending  order  prior  to  summation. 

4)  There  is  considerable  repetition  in  the 
components  from  one  term  to  the  next.  This 
observation  can  be  exploited  to  create  a  fast, 
resource-conserving  algorithm. 
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ABSTRACT:  An  efficient  implementation  of  the 
Generalized  Ray  Expansion  (GRE)  method  for 
computing  the  scattering  of  three-dimensional  (3-D) 
arbitrarily  shaped  deep  cavities  is  studied  in  this 
paper.  Efficiency  is  being  sought  from  two  aspects: 
ray  racing  in  discrete  cavities  and  reflection  from 
individual  patches.  An  improved  algorithm  for 
detecting  intersections  between  a  ray  and  triangular 
patches  has  been  proposed,  which  is  about  2.83 
times  faster  than  the  traditional  algorithm.  Also, 
sectional  algorithm  and  Wavefront  Advancing  and 
Candidate  Narrowing  (WACN)  algorithm  for 
tracing  rays  inside  3-D  cavities  are  proposed  to 
boost  efficiency.  As  to  reflection  from  individual 
patches,  different  local  cavity  reconstruction 
methods  are  being  tested  and  interpolative 
triangular  patches  are  found  to  be  an  efficient 
choice.  Finally,  several  numerical  examples  further 
demonstrate  the  versatility  and  validity  of  our 
approach. 

I.  INTRODUCTION 

Electromagnetic  scattering  from  arbitrarily  shaped 
deep  cavities  is  of  great  importance  in  radar  cross 
section  (RCS)  estimation  of  modem  jet  aircraft  [l]-[8]. 
Because  these  targets  are  usually  composed  of  two 
different  parts,  i.e.,  an  electrically  large,  smooth  varying 
air  duct  and  a  relatively  short,  geometrically  complex 
termination,  methods  suitable  for  one  part  generally 
become  unsuitable  or  even  fail  for  the  other  part.  Due  to 
this  discrepancy,  hybrid  methods  are  often  used  instead 
to  solve  for  different  parts  [8].  In  this  article,  we  shall 
focus  on  efficient  computation  of  the  electrically  large, 
smooth  varying  air  duct.  The  methods  involved 
generally  include  differential  equation-based  methods 
[7j,  integral  equation-based  methods,  waveguide  modal 
analysis  and  high  frequency  asymptotic  methods  [3]. 
Differential  equation-  and  integral  equation-based 
methods  are  accurate  while  much  less  efficient  for  deep 
cavity  problems  due  to  the  prohibitive  amount  of 
memory  and  CPU  requirements.  Furthermore, 
differential  equation-based  methods  suffer  from 
numerical  dispersion  error  for  electrically  large 


problems  and  their  applications  to  deep  cavity 
scattering  are  limited.  Waveguide  modal  analysis  also 
provides  accurate  results  [3],  [4],  but  the  exact 
waveguide  eigenmodes  have  only  been  found  for 
simple  cross  sections.  These  methods  are  most  often 
used  to  give  reference  solutions. 

Because  of  the  smooth  varying  property  of  the  air 
duct  required  by  aerodynamics,  ray  and  beam 
techniques  are  usually  used  for  high  frequency 
asymptotic  methods.  The  early  version  was  the 
Shooting  and  Bouncing  Ray  (SBR)  method  which 
utilizes  Geometric  Optics  (GO)  for  ray  tracing  and 
Aperture  Integration  (AI)  or  Reciprocal  Integral  (RI) 
for  far  field  computations  [l]-[4],  [6].  The  major 
problem  with  the  SBR  method  is  that  it  does  not 
consider  higher  order  effects  -especially  the  field 
diffracted  into  the  cavity  by  the  rim  of  the  open  end. 
Thus  it  generally  provides  an  envelope  but  not  details 
of  the  scattering  pattern.  Gaussian  Beam  (GB)  is 
another  approach  which  instead  traces  Gaussian  beams 
[3],  [5].  Since  the  Gaussian  beam  is  caustic  free  by  its 
nature  and  because  it  considers  fields  diffracted  into  the 
cavity  from  the  open  end,  it  has  much  better  accuracy 
than  the  SBR  method.  But  the  beam  distortion  after  a 
few  reflections  generally  prevents  this  method  from 
deep  cavity  problems.  The  GRE  method  could  be 
thought  of  as  a  combination  of  SBR  and  GB  methods  in 
some  sense  [3].  Based  on  the  sub-aperture  expansion 
techniques  of  the  GB  method,  the  GRE  method  traces 
GO  rays  instead  of  Gaussian  beams  to  improve  the 
beam  distortion  problem.  Since  the  GRE  method 
includes  the  interior  diffraction  by  the  edge  of  the  open 
end,  it  is  also  more  accurate  than  the  SBR  method. 

The  usage  of  the  GRE  method  is  limited  by  the 
necessity  of  tracing  massive  amounts  of  rays.  This  is  in 
turn  related  to  the  modeling  of  a  cavity.  For  simple 
geometry,  analytical  functions  could  be  used  thus  ray 
tracing  is  obviously  not  a  problem.  For  realistic  large 
and  arbitrarily  shaped  cavities,  modeling  a  cavity  with 
ordered  3-D  discrete  points,  which  could  either  be  the 
results  of  physical  measurements  or  generated  by  CAD 
software,  is  of  great  versatility  and  generality. 

When  applying  GRE  methods  to  such  realistic  3-D 
discrete  cavities,  two  essential  issues  need  to  be 
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considered.  The  fist  issue  is  fast  ray  tracing  algorithms 
in  3-D  discrete  cavities.  This  topic  is  rarely  documented 
because  most  3-D  ray  tracing  algorithms  are  designed 
for  3-D  bodies.  We  shall  solve  this  problem  from  two 
aspects:  fastening  the  intersection  test  of  a  ray  and  a 
triangular  patch  and  reducing  the  total  number  of  such 
test  needed.  An  improved  intersection  test  algorithm  is 
proposed  for  the  first  aspect  and  two  other  algorithms,  a 
sectional  algorithm  for  general  cavities  and  a  WACN 
algorithm  for  convex  cavities,  are  proposed  for  the 
second  aspect.  The  second  issue  is  local  cavity 
reconstruction.  There  are  quite  a  few  choices  ranging 
from  simple  triangular  patches  to  complex  Hermitian 
bicubic  patches.  We  shall  study  the  accuracy  of  using 
different  reconstruction  methods  for  reflection 
computation.  We  found  that  an  interpolative  triangular 
patch,  which  is  simple  to  implement  and  highly 
accurate,  was  the  best  choice. 

This  article  is  organized  as  follows:  Section  II  briefly 
introduces  the  GRE  method.  Section  III  discusses  the 
efficiency  issues.  Section  IV  provides  numerical 
examples.  Finally,  some  conclusions  are  drawn  in 
Section  V. 


n.  GRE  METHOD 


In  the  GRE  method,  the  open  end  of  the  cavity  is 
divided  into  multiple  sub-apertures.  The  electric  field 
radiated  by  the  n -th  sub-aperture  is  determined  by  far 
zone  Kirchhoff  approximation  with  the  cavity  wall 
absent.  Cone-shaped  angular  grids  of  ray-tube  are 
launched  from  the  center  of  each  sub-aperture  to 
represent  the  spherical  wave  entering  into  the  cavity. 
By  assuming  a  local  plane  wave  at  the  open  end  and 
using  Physical  Optics  (PO)  to  obtain  the  equivalent 
electric  and  magnetic  currents,  the  electric  field  of  the 
p -th  ray-tube  of  the  n-th  sub-aperture  is  expressed  as 


-A-K * Jf  K xfc xfl,0f)]e'w<S'I 


(1) 

where  k  is  the  wave  number,  fp  =  rp  r  p  and  fn  is  the 
unit  vector  representing  the  direction  of  the  p -th  ray- 
tube,  r' represents  the  location  of  an  equivalent  source 


on  the  sub-aperture,  s'  is  the  unit  surface  normal 


pointing  inwards  the  cavity,  and  E(r)  and  H(r) 
represent  the  incident  electric  and  magnetic  field 
respectively.  The  integration  is  over  the  sub-aperture. 
Note  that  the  far  zone  radiating  field  is  decomposed 


into  the  product  of  a  spherical  wave  and  a  vector  far 
zone  pattern.  Portions  of  the  spherical  wave  could  be 
individually  traced  as  GO  rays.  Incident  field 
information  is  only  contained  in  the  vector  far  zone 

pattern Fn(rnp,Ei).  Ray  tracing  and  the  calculation  of 

Fn(rp ,£()are  independent.  Thus  the  GRE  method 

could  generate  the  result  at  any  incident  angle  in  the 
effective  angular  range  (10°  -  15°  narrower  than  the 
largest  ray  tracing  angle)  with  just  one  ray  tracing. 
Also,  ray  tracing  is  time  consuming  rather  than  memory 
consuming.  The  independency  of  each  ray  tracing 
makes  it  very  suitable  to  utilize  distributed  computer 
systems  because  there  are  virtually  no  communications 
between  different  processes  and  load  balancing  is  easy 
to  handle. 

The  total  transmitted  field  could  be  written  as 

£(')  =  IZ£,f(r/)  (2) 

«=1  p- 1 

To  evaluate Fn ( rp , £,) ,  we  first  express  the  incident 
field  as 


E.<r,‘)  =  P, (3) 

A(r,')  =  AH, (?.')  =  *, xe,(f,')/z0 

where  pr  and  ph  represent  the  directions  of  the  incident 
electric  and  magnetic  field  respectively  and  Z0  is  the 

free  space  wave  impedance.  Establishing  a  local 
coordinate  system,  say  T,  originating  at  the  center  of  a 
sub-aperture  and  in  which  ex,ey  are  any  two  orthogonal 

unit  vectors  tangential  to  the  sub-aperture,  and  in  which 
e.  points  into  the  cavity,  £„(r/,£()  could  be 

A 

decomposed  into  9  and  (p  polarization 


F„  (?/  ,£,.)  =  Fng  (?; ,  £,.  )0n  +  £ ■  (rp ,  £f 


F„eW  >£;)  =  [(p„  +  P„y  cos  6P )  cos  tpp  + 

(Po  -  P,«  cos  6P ) sin  ]I(rp ,  Et ) 


>£,)  =  [(Per  cos 9P  - Phx ) cos <t>p  - 

( Phy +  Pex  cos  O* )  sin  ]/  ( rf ,  Ei ) 


) 

where  6?  and  (p?  represent  the  elevation  and  azimuth 

angles  of  the  axis  of  the  p-th  ray-tube  of  the  n-th  sub¬ 
aperture,  measured  in  the  local  coordinate  system  r\ 
Pc*  ^  Pex  are  components  of  ^in  r'and  similar 


^or  Phx  an<*  Phx  *  exact  form  of  7(r/,£( )  relies  on 
the  shape  of  the  sub-aperture  (see  Ref.  [3]). 

Rays  are  bounced  back  and  forth  inside  the  cavity. 
After  each  reflection,  the  magnitude  is  determined  by 
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(7) 

where  A}  is  the  location  of  the  i-th  reflection, 

■S'  =  |  A}  —  7;_J  I  ,  [DF],_/  is  the  divergence  factor  at  the  (i- 

1  )-th  reflection  location,  [R]  is  the  reflection  matrix  of 
the  cavity  wall  which  could  be  written  in  the  PEC  case 
as 


~E\ 

'-1 

ol 

~E[ 

Er  = 

_L 

Er 

L^//  J 

=  [R]E‘  = 

0 

1_ 

_L 

where  superscripts  r  and  i  denote  the  reflected  and 
incident  wave  respectively,  and  where  and  // 
represent  the  perpendicular  and  parallel  polarization. 
When  a  cavity  is  coated  with  materials,  the  impedance 
boundary  condition  could  be  used  casually  instead  of 
the  PEC  boundary  condition.  Easy  manipulation  of 
boundary  conditions  is  another  advantage  of  ray-based 
techniques  over  other  methods.  The  divergency  factor  is 
determined  by  5  and  the  principal  radii  of  the  curvature 

of  the  wavefront,  say  R\  and  R2i  at  r._x 

[DF]m  =  ,  1 . — •  ,  1  (8) 

Jl  +  s/R,  yjl  +  s/R2 

Note  that  the  reflection  field  is  singular  if  the  i-th 
reflection  is  located  at  the  caustics,  i.e.,  s  =  — /?j 

or  S  =  —  R2 .  The  caustic  problem  is  inherent  to  all  GO- 

based  techniques.  When  it  occurs,  we  chose  to  abandon 
the  ray  being  traced  for  efficiency  considerations. 

Rays  could  exit  from  either  the  front  end  or  the  rear 
end.  In  the  first  case,  the  far  zone  scattering  field  is 
determined  by  the  AI  method.  In  the  second  case,  RI 
could  be  used  to  calculate  the  far  field  contribution  of  a 
ray  tube  directly  without  tracing  it  back.  Without  wall 
losses,  the  cross  section  area  of  the  reflected  ray-cube, 
say  5,  could  be  determined  via  energy  conservation  by 

S0\eo\2  =s|£|2  (9) 

where  S0  is  the  initial  beam  solid  angle  of  the  ray- tube. 

III.  EFFICIENCY  IMPROVEMENTS 

The  major  thrust  of  the  GRE  method  is  to  trace 
massive  amounts  of  ray-tubes  inside  an  arbitrarily 
shaped  cavity.  Usually,  105  -  106  rays  are  expected; 
however,  if  the  axial  length  of  the  cavity  is  about  100>., 
this  number  could  reach  107  -  108.  Therefore,  ray 
tracing  efficiency  is  of  paramount  importance. 

Ray  tracing  is  essentially  a  computer  graphics  topic. 
In  our  context,  its  efficiency  is  determined  rather  by  the 
accuracy  of  scattering  field  computation  than  by  the 
quality  of  graphic  displaying.  We  can  further  divide  the 
ray  tracing  problem  into  two  weakly  coupled  problems: 
determination  of  the  reflection  position  and 


computation  of  the  reflected  field.  In  general,  when  we 
calculate  a  GO  ray  reflection,  we  need  to  reconstruct 
the  local  cavity  from  those  discrete  points  surrounding 
the  reflection  point.  For  a  fixed  set  of  discrete  points,  if 
their  is  a  better  way  to  reconstruct  the  local  surface  so 
that  the  resultant  reflection  calculation  is  more  accurate, 
we  can  use  less  discrete  points  to  model  the  cavity.  The 
efficiency  of  determining  the  reflection  position  in  a 
discrete  cavity  is  essentially  dominated  by  the  number 
of  discrete  points  being  used  to  describe  the  cavity. 
Thus  these  two  problems  are  weakly  coupled  in  this 
sense.  In  general,  we  can  improve  the  overall  efficiency 
by  working  on  each  problem  individually. 

A.  Ray  Tracing 

Ray  tracing  involves  finding  the  reflection  of  a  ray- 
tube.  During  each  ray  tracing,  triangular  patches  were 
used  to  determine  the  reflection  position,  though  in 
some  cases  as  a  preliminary  step.  Basic  ray  tracing 
algorithms  include  two  procedures:  1)  Determination  of 
possible  intersections  of  a  ray  and  all  triangular  patches; 
2)  Sorting  the  distance  between  the  current  position  and 
all  possible  intersections.  The  shortest  distance 
corresponds  to  the  actual  reflection.  We  attempt  to 
improve  the  efficiency  of  each  procedure  in  the 
following. 

A.l  Determine  Intersections 

The  traditional  way  to  determine  a  possible 
intersection  starts  from  calculating  the  intersection  of  a 
ray  and  the  plane  where  a  triangular  patch  is  located 
[9].  Considering  a  ray  originating  at  (x0 ,  y0 ,  z0)  and 
shooting  towards  (kx ,  ky ,  k7 ),  the  ray  function  is  written 
as 

x  =  k](z-z0)  +  x0,y  =  kz(z-z0)  +  y0 

with  kx  —kx  /  kz  and  k2  =  ky  /  kz .  This  requires  two 

multiplications  and  four  summations  (we  do  not 
consider  those  operations  solely  related  to  the  ray 
function  because  they  are  performed  only  once  for  a  ray 
but  not  for  all  triangular  patches  being  tested.).  The 
plane  function  of  a  triangular  patch,  written  as 
ax+by  +  cz  +  d  =  0,  could  be  determined  by  the 
coordinates  of  its  three  vertices  by  solving  a  set  of 
inhomogeneous  linear  equations,  which  requires  36 
multiplications/divisions  and  20  summations.  The 
solution  of  the  intersection  needs  additional  seven 
multiplications/divisions  and  six  summations.  Next,  the 
intersection  is  tested  to  determine  whether  it  is  inside 
the  triangular  patch  or  not.  We  adopt  the  following 
scheme:  1)  Three  vectors  were  constructed  by 
connecting  three  vertices  to  the  intersection.  This 
requires  nine  summations.  2)  All  cross  products 
between  any  two  of  the  three  vectors  were  gathered. 
This  requires  18  multiplications  and  nine  summations. 
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3)  All  dot  products  of  any  two  vectors  obtained  in  step 
two  were  obtained.  This  requires  nine  multiplications 
and  six  summations.  If  the  intersection  is  inside  the 
triangular  patch,  all  cross  products  must  be  in  the  same 
direction  if  being  calculated  with  a  certain  circulation 
order.  Otherwise,  their  must  be  one  cross  product  with  a 
sign  different  from  the  others.  In  fact,  if  we  find  that 
two  dot  products  have  different  signs,  we  can  reach  a 
conclusion  immediately.  This  step  requires  27 
multiplications  and  24  summations  in  the  worst  case. 
Thus  in  the  traditional  method,  totally  72 
multiplications/divisions  and  64  summations  are  needed 
in  the  worst  case  (136  flops). 


Fig.  1.  Comparison  of  traditional  and  proposed  intersection 
determination  approaches. 


The  key  to  efficiency  improvements  is  to  bring  3-D 
operations  to  two-dimensional  (2-D)  operations.  To  do 
so,  we  first  project  the  three  vertices  of  a  triangular 
patch  onto  the  jr-y  plane  of  another  coordinate  system 
originating  at  (%  y0 ,  zo)  and  whose  z- axis  coincides 
with  the  direction  of  the  ray.  This  requires  18 
multiplications  and  21  summations.  Then  all  cross 
products  between  any  two  of  three  2-D  vectors  obtained 
in  the  first  step  are  calculated.  If  the  intersection  is 
inside  the  triangular  patch,  all  three  cross  products  must 
be  of  the  same  sign  when  being  calculated  with  a 
certain  circulation  order.  Otherwise,  their  must  be  one 
cross  product  with  a  sign  different  from  the  others.  This 
step  requires  six  multiplications  and  three  summations. 
Therefore,  the  proposed  scheme  requires  a  total  of  24 
multiplications  and  24  summations  (48  flops). 
Compared  with  the  traditional  approach,  this  algorithm 
needs  35%  less  flops  and  hence  is  2.83  times  faster. 
Moreover,  this  algorithm  is  more  accurate  and  robust 
because  it  does  not  involve  any  division  operation. 

A.2  Ray  Tracing  In  Cavities 

As  has  already  been  pointed  out,  the  heart  of  any  ray 
tracing  algorithm  is  sorting  and  the  key  to  efficiency 
improvements  is  exploiting  data  coherence  [9],  [10].  An 
efficient  algorithm  is  typically  achieved  by  avoiding 
expensive  intersection  computation  as  much  as  possible 
and  by  sorting  the  least  possible  amount  of  intersections 
or  no  such  sorting  at  all.  If  all  3-D  discrete  points  are 


given  without  any  coherence,  it  would  be  hard  to 
improve  the  efficiency. 

Let  us  assume  that  all  discrete  points  are  given  in  m 
consecutive  cuts  along  the  z-axis  and  let  us  call  them  z- 
cut.  Discrete  points  in  one  z-cut  form  a  polygon  (z- 
polygon)  and  those  in  adjacent  cuts  form  a  section  of 
the  cavity  when  connected.  The  whole  cavity  is  formed 
by  (m  -1)  such  sections,  e.g.  Fig.5  and  Fig.6.  With  this 
model,  we  can  search  for  the  next  reflection  section  by 
section  from  where  the  current  reflection  is  and  along 
either  positive  or  negative  z-directions,  depending  on 
the  direction  of  the  ray.  In  this  manner,  the  first 
intersection  must  be  the  actual  reflection  and  no  sorting 
is  needed  at  all.  We  call  this  method  the  sectional 
algorithm. 

In  fact,  3-D  discrete  points  are  either  specified  by 
physical  measurement  or  generated  by  CAD  software. 
It  would  be  natural  to  require  them  to  be  generated  in 
the  above  manner.  For  those  models  which  are  different 
and  can  not  be  regenerated,  we  may  run  a  pre¬ 
processing  program  to  reform  them.  In  the  following, 
we  shall  assume  that  such  a  model  is  always  available. 

The  sectional  algorithm  totally  avoids  sorting,  but  the 
number  of  intersection  computations  could  still  be 
large.  We  can  further  improve  the  performance  for 
convex  cavities  with  the  following  Wavefront 
Advancing  and  Candidate  Narrowing  (WACN) 
algorithm.  This  algorithm  starts  from  computing  the 
intersection  of  a  ray  and  a  z-cut  and  determining 
whether  the  intersection  is  inside  the  z-polygon.  If  the 
intersection  is  out  of  the  current  z-polygon  but  is  inside 
the  previous  z-polygon,  it  must  be  reflected  by  the 
section  formed  by  these  two  z-polygons.  To  test 
whether  an  intersection  is  inside  a  z-polygon,  we  need 
to  specify  a  gauge  point  for  each  z-polygon.  To 
understand  the  role  of  gauge  points,  we  notice  that  any 
2-D  line  (formed  by  adjacent  vertices  of  a  z-polygon) 
divides  a  2-D  space  (or  hyper  plane)  into  two  half¬ 
spaces.  When  all  points  in  one  half-space  are 
substituted  into  the  line  equation,  the  results  must  bear 
the  same  sign  [10].  If  the  line  equation  is  adjusted  such 
that  any  point  from  the  interior  of  the  polygon  yields  a 
positive  (or  negative)  sign  when  being  substituted  into 
the  line  equation,  we  can  determine  whether  an 
intersection  is  inside  a  z-polygon  or  not.  A  gauge  point 
serves  this  purpose  and  it  could  be  any  point  inside  a  z- 
polygon.  A  2-D  line  equation  ax+by+cz  +  d  =  0 
could  be  solved  by  the  coordinates  of  its  two  vertices 
with  six  multiplications/divisions  and  three 
summations.  The  gauging  and  the  determination 
procedures  totally  require  four  multiplications  and  four 
summations.  Thus  for  one  possible  triangular  patch,  ten 
multiplications/divisions  and  seven  summations  (17 
flops)  are  required.  This  is  about  35%  of  the  proposed 
intersection  tests  needed.  Since  the  current  wavefront 
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advances  consecutively,  we  call 
Advancing. 

candidate  points 


it  Wavefront 


candidate  points 


intersection 


previous  z-cut 


current  z-cut 


Fig.  2.  The  local  coordinate  systems  and  candidate  points  that 
the  candidate  triangular  patch  must  contain. 


If  a  reflection  is  about  to  happen  in  a  section,  the 
candidate  triangular  patches  which  are  possible  for 
actual  reflection  could  be  further  narrowed  down.  This 
is  accomplished  through  the  following  steps:  1)  In  the 
previous  z-cut  (where  the  intersection  is  inside  the  z- 
polygon),  construct  a  Cartesian  coordinate  system 
whose  x'-axis  is  the  projection  of  the  ray  on  the  z-cut 
and  whose  y’-axis  is  perpendicular  to  the  x'-axis.  2) 
Transform  all  vertices  of  the  z-polygon  to  this  new 
coordinate  system  and  only  compute  their  y’ 
components.  3)  Check  the  signs  of  all  y*  components 
consecutively.  If  two  adjacent  y*  components  are  of 
opposite  sign,  record  their  indices.  4)  Calculate  the  x * 
components  of  the  two  pairs  of  points  obtained  in  step 
three.  The  candidate  triangular  patch  must  contain  the 
pair  of  points  which  both  have  positive  xy  components. 
5)  Repeat  steps  one  to  three  for  the  current  z-cut  (where 
the  intersection  is  out  of  the  z-polygon).  6)  Calculate 
the  xf  components  of  the  two  pairs  of  points  obtained  in 
step  five.  The  candidate  triangular  patch  must  contain 
the  pair  of  points  closest  to  the  x'  -axis.  After  the  above 
steps,  only  those  triangular  patch  (not  necessarily  two 
patches)  containing  the  two  pairs  of  points  in  the 
current  and  previous  z-cuts  are  possible  for  the  actual 
reflection.  The  candidate  triangular  patches  are 
narrowed  down  and  we  call  this  step  Candidate 
Narrowing.  Note  that  obtaining  either  the  x ' 
components  or  the  y9  components  only  requires  two 
multiplications  and  one  summation,  thus  expensive 
intersection  tests  are  replaced  by  these  simple 
operations. 


B.  Elementary  Reflection 

In  this  subsection,  the  effect  of  using  different  cavity 
reconstruction  methods  for  computing  the  reflection 
will  be  discussed. 


3^-  x  9X  circular  waveguide. 


Fig.  4.  Ocpq,  of  10  7,  x  10  X  circular  waveguide  calculated 
by  using  different  reconstruction  methods. 


B.l  Simple  Triangular  Patch 

In  this  approach,  each  triangular  patch  is  considered 
as  a  simple  plane.  Since  the  principal  radii  are  infinite, 
the  caustic  problem  does  not  exit.  In  general,  more 
triangular  patches  should  be  used  if  better  accuracy  is 
required.  To  study  the  convergency,  we  calculate  the  00 
polarized  mono-static  RCS  of  a  circular  waveguide. 
The  waveguide  is  of  3%  in  diameter  and  9X  in  length. 
Different  numbers  of  triangular  patches  per  section, 
which  are  in  turn  represented  by  the  number  of  discrete 
points  per  circle,  is  being  used.  The  results  are  shown  in 
Fig.3.  We  observe  that  when  only  36  points  per  circle 
(72  patches  per  section)  are  used,  there  are  significant 
errors  for  all  angles.  As  the  number  of  patches  is 
increased,  the  performance  improves  and  180  points  per 
circle  yields  errors  within  3dB/  X2  compared  with  360 
points  per  circle.  Fig.4  shows  the  qxp  polarized  mono¬ 
static  RCS  of  a  circular  waveguide  of  10X  in  diameter 
and  10A,  in  length.  180  points  per  circle  are  used  for 
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reconstruction  by  simple  triangular  patches.  Ray  tracing 
is  confined  within  45°  and  the  effective  angle  is  up  to 
30°  -  35°  according  to  theory  [3].  It  is  observed  that  the 
results  agree  with  those  obtained  by  Modal  analysis 
well  up  to  35°. 

Using  more  triangular  patches  makes  ray  tracing  less 
efficient.  In  the  following,  we  shall  explore  other 
possibilities  with  better  performances. 


B.2  Coons  Patch 

Since  the  inaccuracy  with  simple  triangular  patches  is 
caused  by  the  assumption  of  infinite  principal  radii,  it  is 
natural  to  consider  using  surfaces  with  curvature.  Here 
we  choose  Coons  patch.  The  Coons  patch  belongs  to 
the  family  of  Hermitian  bicubic  parametric  patches.  It 
only  uses  the  information  on  its  four  comers  to 
determine  the  parameters.  To  further  introduce  this 
method,  let  us  denote  a  Coons  patch  with  two 
parameters  (w,  v)  as  R(u,  v)  =  { x(u ,  v),  y(u,  v),  z(m,  v)j, 
with  0  <  u  <  1  and  0  <  v  <  1.  If  uv  is  used  as  an 
abbreviation  for  R(u ,  v),  then  00,  01,  10,  11  represent 
the  four  comers  respectively;  00«,  01«,  10«,  llw,  00», 
01  v,  10v,  1  lv  represent  the  first  order  tangential 
derivatives  at  each  comer;  and  00//v,  01m-,  10«r,  11  uv 
represent  the  second  order  tangential  derivatives  at  each 
comer,  which  are  also  called  twists.  A  Coons  patch  is 
then  expressed  as 

uv  =  (u\u2,u,  1) -[H]  [M]-[H]t  •  (v\ v2, v,l)r 


where 


(10) 


The  great  advantage  of  the  Coons  patch  (as  with 
other  Hermitian  patches)  is  that  when  two  adjacent 
Coons  patches  are  constructed  separately,  the  first  order 
continuity  (C1  continuity)  across  the  patch  edges  is 
guaranteed.  Thus  we  can  construct  a  Coons  patch 
whenever  needed  without  considering  the  global  C1 
continuity.  Compared  with  triangle  patches,  which  are 


of  C°  continuity,  C1  continuity  is  preferred  in  computer 
graphics  because  of  more  realistic  results. 

On  the  other  hand,  the  Coons  patch  is  a  cubic 
function  of  each  of  its  parameter  u  and  v.  This  property 
causes  unnecessary  surface  twists  and  it  is 
disadvantageous  when  being  used  to  calculate  the 
reflection  and  the  divergence  factor  [DF];.  Fig.4  shows 
the  <p<p  polarized  mono-static  RCS  of  a  10X.  x  10}. 
circular  waveguide  calculated  by  using  Coons  patches 
with  36  points  per  circle.  We  observe  that  except  for  the 
main  lobe,  the  results  roughly  deviate  from  the 
reference  values  the  most. 

To  find  the  exact  reflection  position,  we  need  to  solve 
a  set  of  linear  and  non-linear  equations  including  Eq. 
(10)  and  the  ray  function.  If  the  Newton  iterative 
method  is  used,  three  to  four  iterations  should  be 
expected  with  good  initial  guesses  and  appropriate 
accuracy  control.  This  procedure  consumes  more  CPU 
time  than  the  computation  of  the  reflection  field  itself. 

B.3  Interpolate  Triangular  Patch 

In  this  approach,  the  reflection  position  is  determined 
by  treating  a  triangular  patch  as  a  simple  plane.  To 
compute  the  reflection  direction  and  the  divergence 
factor  [DF]„  the  triangular  patch  is  assumed  to  have 
curvature.  Its  first  and  second  order  derivatives  at  the 
reflection  position  are  obtained  by  linear  interpolation 
of  those  at  the  vertices.  Compared  with  the  Coons 
patch,  this  approach  not  only  eliminates  surface  twists, 
but  also  simplifies  the  calculation  of  the  reflection. 

Fig.4  also  depicts  the  qxp  polarized  mono-static  RCS 
of  a  10A,  x  10X  circular  waveguide  calculated  by  using 
interpolative  triangular  patches  with  36  points  per 
circle.  As  can  be  seen,  the  results  are  much  better  than 
those  for  Coons  patches  with  the  same  amount  of  points 
per  circle  and  agree  with  the  reference  values  the  best 
(within  effective  angle  35°).  Compared  with  that  of 
simple  triangular  patches,  the  consideration  of 
curvature  improves  the  accuracy  to  higher  degree  at 
large  angles  (20°  above)  than  at  small  angles.  Bearing 
in  mind  that  the  improvements  are  obtained  by  using 
20%  discrete  points,  as  in  the  case  of  simple  triangular 
patches,  we  consider  this  as  our  best  choice. 

IV.  NUMERICAL  EXAMPLES 

Besides  the  example  shown  in  Fig.4  as  a  verification 
of  our  approach,  we  further  show  some  more 
realistically  shaped  examples  to  demonstrate  the 
versatility.  The  first  example  is  a  PEC  cavity  with  a 
slanted  front  end.  The  cavity  is  formed  by  six  sections. 
The  first  section  is  a  slanted  aperture  of  10>.  in  length. 
The  angle  between  the  normal  to  the  aperture,  i.e.  fi , 
and  the  z-axis,  is  45°.  The  angle  between  the  plane,  in 
which  h  and  the  z-axis  are  located,  and  the  x-axis  is  0°. 
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The  second  section  is  10X  in  length  and  its  cross  section 
is  a  square  with  side  lengths  of  10X.  The  third  section  is 
a  6.7X  transition  region  where  the  cross  section  changes 
from  a  square  with  a  side  length  of  1  OX  to  a  circle  of 
10X  in  diameter.  The  fourth  section  is  8.3X  in  length  and 
its  cross  section  is  a  circle  of  1  OX  in  diameter.  The  fifth 
section  is  another  transition  region  of  8.3X  in  length  and 
its  cross  section  changes  from  a  circle  of  10X  in 
diameter  to  a  circle  of  8X  in  diameter.  The  final  section 
is  of  6.7X  in  length  and  its  cross  section  is  a  circle  of  8X 
in  diameter.  The  geometry  is  shown  in  Fig.5  with 
adjusted  axial  ratios.  The  side  view  is  shown  in  Fig.7 
with  real  axial  ratios.  All  figures  are  in  the  unit  of 
wavelength. 


Fig.  5.  Cavity  with  slanted  front  end. 


Fig.  6.  Cavity  with  axial  lofting. 


For  comparison,  we  built  another  model  with  a 
normal  front  end.  The  only  difference  is  that  the  first 
section  has  the  same  cross  section  as  the  second  section, 
i.e.  the  front  end  is  perpendicular  to  the  z-axis.  Both 
cavities  are  terminated  with  simple  PEC  plates. 

We  use  interpolative  triangular  patches  to  reconstruct 
the  cavity  with  72  points  per  cross  section  and  WACN 
algorithm  for  ray  tracing.  Fig.9  and  Fig.  10  show  that 
both  00  and  99  polarized  RCS  at  9  =  0°.  0  takes  a 


positive  sign  when  an  observation  has  a  positive  x 
coordinate  and  negative  sign  otherwise.  As  we  see,  the 
main  lobes  of  a9G  and  ow  of  the  cavity  with  normal 
front  end  do  not  occur  at  the  normal  incidence  but  at 
some  larger  angles.  When  the  front  end  becomes 
slanted,  the  main  lobe  is  close  to  the  normal  incidence 
but  shifts  slightly  towards  negative  0  direction. 


x-axis 

Fig.  7.  Side  view  of  the  cavity  with  slanted  front  end. 

Note  that  the  results  for  the  cavity  with  normal  front 
end  are  not  exactly  symmetrical.  This  is  due  to  the  low 
grid  density  (eight  points  per  X)  being  used  in  aperture 
integration.  The  results  converge  slowly  to  symmetrical 
forms  if  grid  density  becomes  denser.  Without 
exception,  all  aperture  integration  in  this  section  will  be 
performed  with  the  above  grid  density. 

The  second  example  is  a  concave  cavity  with  axial 
lofting  as  shown  in  Fig.6  with  adjusted  axial  ratios.  A 
side  view  with  real  axial  ratios  is  depicted  in  Fig.  8.  The 
axis  is  described  by  the  following  function  with  z  as  a 
parameter 

x  =  0,  y  =  2(X  -  cos (nz  /  10(U)) 

Each  cross  section  is  formed  by  two  parts.  The  shorter 
one  is  an  arc  of  a  circle  of  5X  in  radius.  The  longer  one 

is  described  by  a  curve  5  A(^3A  -  sin 2  yr  +  cos  iff)  with 

47i/3  <  9  <  871/3.  Both  are  centered  at  the  axis.  The 
cavity  is  also  terminated  by  a  PEC  plate. 

We  use  interpolative  triangular  patches  to  reconstruct 
the  cavity  with  24  points  per  arc,  and  sectional 
algorithm  for  ray  tracing.  Fig.  11  and  Fig.  12  reveal  that 
the  RCS  of  00  and  99  polarization  at  9  =  0°  and  9  = 
90°.  At  9  =  0°,  0  takes  a  positive  sign  when  an 
observation  has  a  positive  jc  coordinate  and  a  negative 


WANG,  LI,  WANG,  ZHU:  GRE  TECHNIQUES  FOR  3D  CAVITIES 


30 


sign  otherwise.  A  similar  convention  holds  at  <p  =  90°. 
The  RCS  at  <p  =  0°  has  main  lobes  at  normal  incidence. 
The  weak  asymmetry  is  also  caused  by  insufficient  grid 
density  in  aperture  integration.  At  <p  =  90°,  the  main 
lobes  of  both  polarization  shift  toward  the  negative 
direction.  This  corresponds  to  the  direction  where  the 
termination  could  be  illuminated  directly.  Note  that 
there  are  actually  two  main  lobes  for  the  80 
polarization  and  its  scattering  is  much  stronger  than  that 
of  the  cp(p  polarization. 


Fig.  8.  Side  view  of  the  cavity  with  axial  lofting. 
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Fig.  10.  Comparison  of  at  (p  =  0°  of  cavities  with  slanted 
and  normal  aperture. 


Fig.  1 1.  Coo  and  at  cp  =  0°  of  the  cavity  with  axial  lofting. 
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V.  CONCLUSION 

In  this  article,  we  have  discussed  several  efficiency 
considerations  of  the  GRE  method  for  computing 
electromagnetic  scattering  from  3-D  arbitrarily  shaped 
deep  cavities.  An  improved  algorithm  for  testing  the 
intersection  of  a  ray  and  a  triangular  patch  is  proposed, 
which  is  2.83  times  faster  than  the  traditional  approach. 
Two  efficient  algorithms  for  ray  tracing  in  3-D  discrete 
cavities  -  the  sectional  algorithm  and  the  WACN 
algorithm  -  are  also  proposed.  The  WACN  algorithm 
further  boosts  the  efficiency  by  2.83  times  for  convex 
cavities.  The  effects  of  using  different  reconstruction 
methods  are  explored.  Numerical  examples  further 
show  the  validity  and  versatility  of  our  approaches. 
Future  work  should  address  the  implementation  of  these 
approaches  on  distributed  computer  systems. 
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Abstract: 

We  discuss  and  demonstrate  by  measurements  and 
computations  the  relation  between  electromagnetic 
bandgap  surfaces  (EBG)  used  to  realize  artificial 
magnetic  conductors  and  the  so-called  soft  and  hard 
surfaces  in  electromagnetics,  with  respect  to  their  STOP 
and  GO  characteristics  for  surface  waves.  We  show 
how  the  main  characteristics  of  such  surfaces  can  be 
modeled  by  using  ideal  surfaces  representing  perfect 
magnetic  conductors  (PMC)  and  PEC/PMC  strip  grids. 
Unfortunately,  commercial  codes  do  not  allow  such 
modeling  for  general  shapes  of  the  surfaces. 

1.  Introduction  to  EBGs,  AMCs  and  soft 
and  hard  surfaces 

In  the  last  few  years,  there  has  been  much  research  on 
using  periodic  structures  to  make  new  materials  for 
application  in  the  design  of  antennas  and  microwave 
components.  Such  materials  are  most  often  referred  to 
as  photonic  bandgap  structures  (PBG),  electromagnetic 
bandgap  structures  (EBG),  or  electromagnetic  crystals. 
In  some  of  this  work,  the  main  concern  is  the  surface 
characteristics  of  the  periodic  structures,  in  the  sense 
that  the  goal  is  to  obtain  a  high  surface  impedance  [1], 
or  to  remove  surface  waves  from  a  dielectric  substrate. 
This  latter  characteristic  is  herein  referred  to  as  a 
STOP-characteristic.  A  high  impedance  surface 
represents  an  artificial  magnetic  conductor  (AMC). 
AMCs  have  the  characteristic  that  electric  current 
sources  such  as  dipoles  can  be  located  at  the  surface 
and  still  radiate  well.  Thereby  very  low  profile  antennas 
can  be  realized  [2].  Other  applications  of  AMCs  are  to 
realize  waveguides  that  can  support  TEM  waves  [3], 
which  in  this  paper  is  referred  to  as  a  GO-characteristic. 
Actually,  the  waves  should  be  regarded  as  quasi-TEM 


waves,  as  they  only  appear  at  the  frequency  where  the 
surface  impedance  is  infinite. 

A  classical  relative  to  the  above-mentioned  surfaces  is 
the  transversely  corrugated  surface.  These  types  of 
surfaces  were  originally  been  used  as  chokes  to  reduce 
coupling,  see  e.g.  the  recent  papers  [4]  and  [5],  and  in 
corrugated  horn  antennas  that  have  found  so  many 
applications  in  large  reflector  antennas.  In  1987  the 
relation  between  the  corrugated  surfaces  and  the  so- 
called  soft  and  hard  surface  described  in  diffraction 
theory  and  acoustics  was  discovered  (in  acoustics  the 
soft  and  hard  surfaces  are  actually  soft  and  hard  to 
touch).  In  terms  of  the  boundary  conditions  of  the 
fields,  an  edge  in  a  transversely  corrugated  surface 
should  be  analyzed  as  an  edge  in  a  soft  surface  by  using 
the  soft  diffraction  coefficient,  both  in  E-plane  and  Pi- 
plane  (STOP-characteristics  in  both  planes).  In 
comparison,  for  a  normal  smooth  conductor  the  edge  is 
soft  in  H-plane  and  hard  in  E-plane.  This  fact  was  used 
to  define  soft  and  hard  surfaces  in  electromagnetics  [6]- 
[7]  in  terms  of  surface  impedances  and  the  boundary 
conditions  in  E-  and  H-  planes.  In  short,  the  soft  and 
hard  surfaces  are  anisotropic.  The  soft  surface  behaves 
like  a  perfect  electric  conductor  (PEC)  in  H-plane  and 
as  a  perfect  magnetic  conductor  (PMC)  in  E-plane,  and 
visa  versa  for  the  hard  surface,  see  Figure  I.  The 
preferred  illustration  of  an  ideal  soft-hard  surface  today 
is  a  PEC/PMC  strip  grid  with  zero  strip  period,  or  in 
other  words,  a  surface  with  electric  and  magnetic 
conductivity  in  one  (and  the  same)  direction  only,  see 
Figure  2a.  Note  that  the  surfaces  in  the  figure  are  color- 
coded  with  blue  meaning  PEC  and  gold  meaning  PMC. 
The  PEC/PMC  strip  grid  represents  a  hard  surface  when 
the  strips  are  oriented  in  the  same  direction  as  the  wave 
propagates  (longitudinal  strips),  and  a  soft  surface  when 
they  are  oriented  orthogonal  to  this  direction  (transverse 
strips).  The  concept  of  soft  and  hard  surfaces  has  later 
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been  generalized  to  arbitrary  orthogonal  PEC/PMC 
anisotropy  [8]-[9],  but  this  will  not  be  treated  here. 


Surface  Name 

new  classical 

Polari; 

VER 

nation 

HOR 

mm 

GO 

STOP 

IHH  PMC 

STOP 

If  IIP  II  ffl 

STOP 

STOP 

Mil  ;!  II ,  m , 

|p$pj|a*na£ 

GO 

vsjr  -i J! 

GO  surfaces:  Enhances  propagation  of  waves  along. 
STOP  surface:  Stops  propagation  of  waves  along. 

PEC  =  Perfect  Electric  Conductor 
PMC  =  Perfect  Magnetic  Conductor 

AMC  =  Artificial  Magnetic  Conductor 
PBG  =  Photonic  Bandgap  Material 
EBG  =  Electromagnetic  Bandgap  Material 
EMXtals  =  Electromagnetic  Crystals. 

Figure  1.  Characteristics  of  different  types  of  surfaces 
with  respect  to  propagation  of  surface  waves  for 
different  polarizations  (top),  and  explanation  of 
abbreviations  (bottom).  Note  that  we  here  use  the  term 
surface  waves  in  an  extended  sense,  so  that  they  even 
include  grazing  waves  along  a  PEC  surface  (behaving 
like  guided  surface  waves  at  cut-off). 


Figure  2b.  Three  modern  realizations  in  terms  of  metal 
strips  on  grounded  dielectric  substrates.  The  upper  case 
consists  of  close  and  narrow  metal  strips,  in  which  case 
the  thickness  of  the  dielectric  layer  determines  the  fre¬ 
quency  at  which  the  surface  is  soft  or  hard.  The  lower 
case  consists  of  wide  metal  strips,  in  which  case  the 
strip  width  determines  the  soft/hard  frequency.  The 
performance  is  considerably  improved  if  the  strip  is 
grounded  with  metal  posts  or  vias,  as  shown.  In  the 
middle  case  the  metal  posts  or  vias  are  located  along 
one  edge  of  the  strips,  in  order  to  reduce  the  width  to 
half.  The  effective  permittivity  in  the  formula  is  for  the 
transverse  strips  (soft  case)  given  by  £eff  =  £r  for  the 
lower  two  geometries.  For  the  hard  case  and  the  upper 
soft  case  they  are  given  by  Seff  =  Sr  - 1 . 

2.  EBG  realizations  of  soft  and  hard 
surfaces 


Figure  2a.  PEC/PMC  strip  representations  of  ideal  soft 
and  hard  surfaces.  The  red  and  green  wave-shaped 
arrows  represent  the  direction  of  propagation  of  the 
waves  that  makes  the  PEC/PMC  strip  surface  soft  and 
hard,  respectively. 


The  original  realizations  of  the  soft  and  hard  surfaces, 
which  were  initially  proposed  and  studied,  were 
transverse  corrugations  (soft),  dielectric-filled 
longitudinal  corrugations  (hard),  and  dielectric 
substrates  with  transverse  (soft)  or  longitudinal  (hard) 
metal  strips.  These  realizations  give  a  thick  surface 
compared  to  common  EBGs,  because  the 

thickness  h  =  Xj  A^S^ ,  where  £eff  is  defined  in  the 

caption  of  Figure  2b.  However,  new  and  much  thinner 
EBG-inspired  realizations  are  readily  available  as 
shown  in  Figure  3,  even  if  the  detailed  performance  has 
not  yet  been  investigated.  The  AMCs  work  also  without 
metal  posts  ,  but  the  performance  is  reduced.  In  the 
same  way  the  strip-loaded  soft  and  hard  surfaces  work 
without  metal  posts  or  vias,  but  the  performance  is 
much  better  with  them.  The  strip-loaded  soft  and  hard 
surfaces  suffer  from  severe  problems  with  strip  modes 
that  are  effectively  killed  with  close  metal  connections 
to  the  ground  such  as  posts  or  vias. 
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3.  Applications  of  STOP  Surfaces 

An  important  application  of  EBGs  and  AMCs  today  is 
to  reduce  the  sidelobes  in  E-plane  (STOP- 
characteristic)  for  aperture  [12]  and  microstrip  antennas 
[11].  The  principle  of  operation  is  readily  explained  in 
terms  of  the  table  in  Figure  1.  From  this  table  it  is  also 
clear  that  a  soft  surface  will  provide  sidelobe  reduction 
(STOP-characteristic)  of  antennas  in  any  plane  and  for 
any  polarization,  such  as  the  cases  treated  in  [13]. 
Different  realizations  of  soft  surfaces  for  sidelobe 
reduction  are  studied  in  [14].  Recent  papers  also 
propose  metal  strips  in  combination  with  an  AMC 
(making  it  a  soft  surface)  to  reduce  coupling  [15]. 


4.  Application  of  GO  surfaces  (quasi-TEM 
waveguide  and  hard  horn) 

Another  application  of  AMCs  and  EBGs  is,  as  already 
mentioned,  TEM  waveguides  [3].  By  using  a  PMC  on 
the  E-plane  walls,  a  rectangular  waveguide  can  support 
a  TEM  wave.  Actually,  if  all  the  waveguide  walls  are 
made  of  a  hard  surface  the  TEM  performance  is  better, 
because  a  TEM  wave  of  any  polarization  can  propagate 
in  such  a  hard  waveguide  of  arbitrary  cross  section. 
This  was  described  already  in  the  original  papers  on 
soft  and  hard  surfaces.  The  first  simulation  of  the  field 
solution  in  a  hard  quasi-TEM  waveguide  was  done  by  a 
FDTD  code  in  [16].  The  simulations  were  done  both  for 
realizations  in  terms  of  corrugations  and  strips,  and  the 
latter  case  included  also  a  study  of  an  homogenized 
asymptotic  model  for  the  strip  grid  [17],  Two  simulated 
cases  are  shown  in  Figure  2.  We  see  the  uniform  field 
distribution  over  the  cross  section  and  evanescent  fringe 
fields  around  the  strips  (Figure  2a).  In  practice,  the 
TEM  waves  of  these  ideal  guides  can  only  be  present  at 
the  center  frequency  where  the  surfaces  have  close  to 
ideal  performance.  Still,  their  bandwidth  is  sufficient  to 
be  attractive  for  use  in  cluster-fed  multi-beam  antennas 
for  Ka  band  multimedia  satellites.  They  are  also  studied 
at  several  places  for  use  in  quasi-optical  grid  amplifiers 
in  the  millimeter  wave  region  (see  e.g.  [21]).  Some 
attempts  to  realize  hard  homs  with  uniform  aperture 
distribution  were  made  more  than  10  years  ago,  but  at 
that  time  the  numerical  techniques  were  not  developed 
sufficiently  to  control  of  the  performance.  Today,  it  is 
possible  to  analyze  hard  homs  very  accurately  with 
commercial  software  based  on  FEM  or  FDTD 
approaches.  The  homs  can  be  analyzed  more  time- 
efificiently  by  the  special  mode  matching  technique 
used  in  [22]-[23].  This  makes  use  of  asymptotic  strip 
and  corrugation  boundary  conditions  [17]  to  simplify 
the  analytical  modal  expansion  in  each  cylindrical 
section  (see  Figure  4). 


Hard  surfaces  can  also  be  used  to  let  waves  GO  past  or 
through  obstructions  without  generating  blockage  [24], 
and  to  let  them  GO  along  a  metal  cylinder  [25]. 


Asymptotic  strip  boundary 
condition 


Dsttfibutiw  of  Hy-ficld 


Figure  3a.  Geometry  (upper)  and  computed  H-fleld 
distribution  (lower)  in  dual-polarized  quasi-TEM  hard 
waveguide  realized  by  strip-loaded  dielectric  substrate. 
This  case  has  been  computed  by  using  the  asymptotic 
strip  boundary  condition  [17], 


Strip-loaded  dielectric  walls 
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Figure  3  b.  Geometry  (upper)  and  computed  H-field 
distribution  (lower)  in  dual-polarized  quasi-TEM  hard 
waveguide  realized  by  modeling  each  metal  strip  of 
finite  thickness .  The  computations  have  been  done  by 
FDTD  as  explained  in  [16]. 


35 


ACES  JOURNAL,  VOL.  18,  NO.  1,  MARCH  2003 


Figure  4 .  Hard  horn  geometry  approximated  as 
cascaded  sections  of  cylindrical  corrugated  sections , 
for  analysis  by  mode  matching  and  generalized 
scattering  matrices  [23].  The  fields  in  each  corrugated 
cylindrical  section  is  expanded  in  modes  by  making  use 
of  the  asymptotic  corrugation  boundary  condition  [17]. 
A  similar  mode  matching  approach  has  been  developed 
for  strip-loaded  hard  horns. 


5.  About  EM  modelling  of  complex  surfaces 

The  needs  for  future  research  on  EM  modelling  in 
connection  with  AMCs  and  soft  and  hard  surfaces  can 
be  divided  in  four  categories: 


They  have  also  been  implemented  in  software  based  on 
the  moment  method  for  cylindrical  structures  of 
arbitrary  cross  section,  such  as  e.g.  [30].  Commercial 
codes  can  normally  not  make  use  of  homogenized 
boundary  conditions.  The  impedance  boundary 
condition  is  also  a  homogenized  boundary  condition 
that  is  applied  to  such  surfaces  [31].  However,  it  is  not 
very  convenient  when  modelling  soft  and  hard  surfaces 
because  it  will  vary  with  angle  of  incidence.  The 
asymptotic  boundary  conditions  are  much  more 
accurate. 

C.  Studies  based  on  exact  modeling. 

The  exact  modeling  of  the  surfaces  in  all  details  is 
always  desirable,  but  it  can  rarely  be  done  due  to  the 
computational  effort  involved.  The  segment  size  often 
becomes  too  small  to  accurately  model  periodicity  and 
thereby  the  computation  time  and  memory  becomes 
larger  than  for  smooth  surfaces. 

D.  Ray  techniques. 

The  computational  effort  with  FEM,  FDTD  and 
moment  method  for  large  structures  can  be  reduced  by 
using  ray  techniques.  These  will  not  be  mentioned  here, 
except  referring  to  one  of  the  papers  in  this  area  [32]. 


A.  Studies  based  on  ideal  surface  models. 

Existing  literature  does  not  contain  any  theoretical 
studies  of  characteristics  of  ideal  PMCs  or  soft  and  hard 
surfaces,  such  as  e.g.  waveguide  solutions  and  Green’s 
functions.  This  would  be  desirable  in  order  to  foresee 
possible  applications.  Some  Green’s  functions  are 
available  in  [26]  focusing  mainly  on  studies  of  surface 
waves.  Commercial  codes  can  often  model  infinite 
plane  PMCs,  but  finite  PMCs  and  PEC/PMC  strip 
models  are  not  included.  Sometimes  the  latter  can  be 
modelled  approximately  by  locating  PEC  strips  on  an 
infinite  AMC  surface. 

B.  Studies  based  on  homogenized  boundary  cond. 

The  major  performance  of  a  surface  in  different 
applications  can  most  effectively  be  found  if  the 
periodicity  of  the  surface  can  be  removed  by 
homogenization  of  the  boundary  condition.  The 
asymptotic  strip  and  corrugation  boundary  conditions  in 
[17]  are  examples  of  homogenized  boundary 
conditions.  They  have  already  been  used  in  [16],  [18]- 
[20],  [23]  and  [26].  The  homogenized  strip  and 
corrugation  boundary  conditions  can  also  be  used  to 
derive  analytical  solutions  of  realized  soft  and  hard 
waveguides,  such  as  the  circular  strip-loaded  guide  in 
[27].  They  have  been  used  in  FEM  software  [28],  and 
they  have  been  implemented  in  algorithms  and  software 
based  on  plane  wave,  cylindrical  mode  and  spherical 
mode  analysis  of,  respectively,  planar,  circular 
cylindrical  and  spherical  multilayer  structures  [29]. 


6.  Experimental  illustration  of  relation 
between  AMC-type  EBGs  and  soft  and 
hard  surfaces 

In  this  section,  we  will  illustrate  the  relation  between 
AMCs  and  soft  and  hard  surfaces  by  taking  an  AMC 
and  transforming  it  to  a  soft/hard  surface  by  using  metal 
tape  in  the  form  of  narrow  strips,  see  Figure  5.  The 
EBG  surface  was  obtained  from  Ericsson  AB.  Some 
previous  work  on  it  was  presented  in  [33].  It  consists  of 
two  patch  layers  of  Sievenpiper  mushroom  type  [1]  for 
900  MHz  applications.  All  the  patches  on  both  layers 
are  connected  to  the  metal  ground  plane  by  metal  posts. 
We  choose  to  measure  coupling  for  both  horizontal 
(HOR)  and  vertical  (VER)  polarizations  between  two 
quarterwave  vertical  monopoles  and  two  horizontal 
parallel  halfwave  dipoles,  respectively.  The  results  in 
Figures  5c  and  5d  show  the  coupling  after  we  have 
corrected  for  the  mismatch  of  the  two  antennas  by  using 
the  following  formula: 


C„=101og 


(l-|S„|2)(l-|S22|2) 


The  formula  represents  the  coupling  that  would  be 
measured  if  we  match  both  antennas  and  neglect 
multiple  interactions  between  them  (i.e.  coupling  to  the 
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Figure  5a.  Photo  of  the  dual-layer  mushroom-type 
AMC  of  Sievenpiper  type. 


Frequency  (MHz) 


Figure  5c.  Measured  net  transmissions  between  two 
vertical  quarterwave  monopoles  on  4  different  surfaces: 
The  metal  plate,  the  original  AMC,  the  strip-loaded 
hard  AMC  (shown  in  Figure  5b),  and  the  strip-loaded 
soft  AMC. 


Figure  5d.  Measured  net  transmission  between  two 
horizontal  halfwave  dipoles  located  9mm  above  the  4 
different  surfaces. 


Figure  5b.  Two  monopoles  on  the  dual-layer 
mushroom-type  AMC  in  Figure  5a  that  has  been 
transformed  to  a  hard  surface  by  providing  it  with 
longitudinal  metal  strips  made  from  Aluminum  tape.  If 
the  strips  are  taped  to  the  surface  in  the  opposite 
transverse  direction ,  the  surface  becomes  soft. 


neighboring  dipole  and  back  again).  This  is  satisfied  if 
the  distance  between  them  is  sufficiently  large.  We 
refer  to  (1)  as  the  net  coupling.  The  two  antennas  were 
mounted  at  a  fixed  distance,  and  four  different  cases 
were  measured:  metal  ground  plane  of  same  size  as  the 
EBG  (corresponding  to  a  PEC),  the  original  EBG 
surface  working  as  an  AMC,  the  EBG  with  longitudinal 
metal  strips  on  it  (corresponding  to  a  hard  surface,  see 
Figure  5b),  and  the  EBG  with  transverse  metal  strips  on 
it  (corresponding  to  a  soft  surface).  The  measured 
results  show  that  both  the  AMC  and  the  soft  surface 
have  a  clear  bandgap  (i.e.  STOP  band)  for  VER 
polarization.  For  HOR  polarization,  there  is  a  STOP 
band  below  900  MHz,  but  a  GO  band  above.  Thus,  the 
EBG  surface  acts  as  an  AMC  around  900  MHz  for  VER 
polarization,  and  it  transforms  gradually  into  STOP  one 
at  that  frequency  for  HOR  polarization.  The  PEC 
STOPs  the  waves  for  HOR  polarization  and  lets  them 
GO  for  VER  polarization.  The  soft  surface  effectively 
stops  the  waves  for  both  VER  and  HOR  polarization, 
whereas  the  hard  surface  lets  them  GO.  The  STOP 
characteristics  of  the  soft  surface  are  not  as  good  for 
VER  polarization  as  for  the  AMC  case,  but  this  may  be 
better  if  actual  anisotropic  soft  EBGs  of  the  kind  shown 
in  Figure  2b  are  realized.  It  seems  to  be  possible  to  use 
the  original  EBG  as  a  STOP  surface  for  both  HOR  and 
VER  polarization  slightly  below  900  MHz,  due  to  some 
frequency  shift  between  the  AMC  performances  for  the 
two  polarizations.  The  GO  characteristics  of  the  hard 
surface  are  very  clearly  present  for  both  polarizations. 
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1.  Numerical  results  for  ideal  surfaces 
(PEC,  PMC,  soft  and  hard  surfaces) 

We  have  also  computed  the  net  coupling  between  VER 
monopoles  and  HOR  dipoles  located  on  the  different 
corresponding  ideal  surfaces.  For  this  we  have  used 
three  different  moment  method  codes,  one  commercial 
code  WIPL-D  [34]  and  two  in-house  codes.  The  latter 
are  WireMoM  [35]  and  G2DMULT  [30]. 

1 .  The  commercial  code  models  the  actual  finite  metal 
surfaces  as  an  ideal  PEC,  and  the  currents  on  both  the 
ground  plane  and  the  wires  are  solved  for  using  the 
entire  domain  basis  functions.  The  PMC  modeling 
assumes  an  infinite  ground  plane.  Soft  and  hard 
surfaces  were  modeled  by  locating  metal  strips  on  the 
PMC. 

2.  The  WireMoM  program  calculates  the  current 
distribution  on  any  wire  antenna  by  using  subsectional 
basis  functions.  TTie  PEC  and  PMC  ground  planes  are 
included  by  imaging. 

3.  G2DMULT  uses  a  spectrum  of  two-dimensional 
solutions  solved  by  the  moment  method  to  determine 
the  coupling  between  antennas  with  a  given  cosine¬ 
shaped  current  distribution.  It  can  handle  all  kind  of 
ground  planes:  PEC,  PMC,  soft  and  hard.  The  finite 
widths  of  the  different  ground  planes  are  included, 
which  are  assumed  infinitely  long. 

The  measured  and  computed  results  are  shown  in  the 
Tables  1  and  2.  Some  values  are  missing  in  the  tables 
because  the  codes  could  not  be  used,  or  the  results  were 
unreasonable.  In  particular,  it  was  difficult  to  calculate 
horizontal  dipoles  located  close  to  a  PEC  and  a  soft 
surface,  and  vertical  monopoles  on  a  PMC.  The  reasons 
are  that  the  PEC  and  soft  surface  short-circuits  HOR 
electric  currents  and  make  their  impedances  close  to 
zero,  and  the  impedances  approach  infinity  for  vertical 
monopoles  on  a  PMC.  The  presented  computed  values 
correspond  well  to  the  measurements  although  we 
cannot  of  course  predict  the  frequency  variations  due  to 
the  actual  surfaces  with  such  ideal  models,  but  the 
major  STOP  and  GO  characteristics  are  well  predicted. 


8.  Conclusion 

We  have  demonstrated  the  relation  between  PBG 
surfaces  of  AMC-type  and  the  soft  and  hard  surfaces 
regarding  their  STOP  and  GO  characteristics  with 
respect  to  surface  waves.  We  have  also  shown  that  there 
is  a  large  need  for  being  able  to  model  finite  and 
arbitrarily  shaped  ideal  PMCs  and  soft  and  hard 


surfaces  by  commercial  codes,  but  first  the  numerical 
techniques  must  be  further  developed  to  enable  such 
code  extensions. 

Table  1.  Measured  net  coupling  levels  in  dB  between 
wire  antennas  on  different  surfaces  (ground  planes)  for 
an  antenna  spacing  of  1.26^  in  a  frequency  band  of 
about  50  MHz  around  900  MHz.  VER  means  vertical, 
HOR  means  horizontal,  h  is  height  over  ground.  The 
STOP  cases  with  low  coupling  are  marked  with  bold. 
(Soft  and  hard  are  made  by  strips  on  AMC). 


VER  monopoles 

HOR  dipoles 

HOR  dipoles 

E-plane 

H-plane,  h=0 

H-plane,  h=A,/4 

Metal 

-21 

-28  to  -40 

-22  to  -35 

AMC 

-43  to  -26 

-30  to -10 

-20  to  -1 1 

Soft 

-32  to  -23 

-32  to  -45 

-22  to  -35 

Hard 

-5  to  -14 

-10  to -17 

-11  to -15 

Table  2.  Computed  net  coupling  levels  in  dB  between 
wire  antennas  on  different  infinite  ideal  surfaces  for  an 
antenna  spacing  of  1.26A,  at  900  MHz.  The  different 
results  are  obtained  by  using  different  codes  based  on 
the  moment  method.  VER  means  vertical,  HOR  means 
horizontal,  h  is  the  height  over  ground  and  x  represents 
missing  data.  The  STOP  cases  with  low  coupling  are 
marked  with  bold.  (Soft  and  hard  are  made  by 
PEC/PMC  strips). 


VER  monopoles 

HOR  dipoles 

HOR  dipoles 

E-plane 

H-plane,  h=0 

H-plane,  h=A./4 

PEC 

-22,  -20 

X,  -30,X 

-26,  -26,  -26 

PMC 

* 

1 

<*> 

o 

-20,-20,-21 

-14,-13,-14 

Soft 

-29,  X 

X,  -30,  X 

-26,  -23,  X 

Hard 

-10,X 

-10,X,X 

-13,  X,  X 

Codes 

Used 

G2DMULT, 
Wire  Mom 

G2DMULT, 
WIPL-D, 
Wire  Mom 

G2DMULT, 
WIPL-D, 
Wire  Mom 
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Theoretical  and  Experimental  Investigations  of  the 
Surge  Response  of  a  Vertical  Conductor 

Md.  Osman  Goni,  Hideomi  Takahashi 


Abstract —  This  paper  describes  the  theoretical,  simula¬ 
tion  and  experimental  investigations  of  the  surge  response 
of  a  vertical  conductor,  including  the  effects  of  ground  sur¬ 
face  and  without  ground  surface.  One  of  the  authors  of 
this  paper  derived  the  formula  of  the  surge  impedance  in 
case  with  ground  surface  and  without  ground  surface.  In 
this  research,  these  theoretical  formulas  are  examined  by 
the  simulation  analysis  of  the  vertical  conductor  using  the 
Numerical  Electromagnetic  Code  (NEC-2)  as  well  as  the 
experimental  basis.  The  arrangement  of  the  current  lead 
wire  in  the  vertical  conductor  model  to  be  analyzed  here 
is  verified  with  the  simulation  result  of  the  equivalent  cir¬ 
cuit  model  by  the  Electromagnetic  Transients  Program 
(EMTP). 

Keywords — EMTP,  Lightning  surge,  Numerical  electro¬ 
magnetic  field  analysis,  Tower  surge  impedance,  Vertical 
conductor. 

I.  INTRODUCTION 

REDICTION  of  lightning  surges  is  very  impor¬ 
tant  for  the  design  of  electric  power  systems  and 
telecommunication  systems.  In  particular,  tower  surge 
impedance  is  an  important  factor  in  analysis  of  the  light¬ 
ning  performance  of  transmission  lines.  Therefore,  a 
number  of  experimental  and  theoretical  studies  on  tower 
surge  impedance  have  been  carried  out  [1]— [11] . 

The  first  theoretical  formulation  of  tower  surge 
impedance  was  proposed  by  Jordan  [1].  The  tower  was 
approximated  as  a  vertical  cylinder  having  a  height  equal 
to  that  of  the  tower,  and  a  radius  equal  to  the  mean 
equivalent  radius  of  the  tower.  Theoretical  formulations 
of  tower  surge  impedance  based  on  the  electromagnetic 
field  theory  were  proposed  by  Lund  holm  et  al.  [2],  Wag¬ 
ner  and  Hileman  [3],  Sargent  and  Darveniza  [4]  and  Oku- 
mura  and  Kijima  [5],  considering  effects  of  the  vector  po¬ 
tential  generated  by  the  injection  current  into  the  tower 
only. 

Another  experimental  value  for  actual  transmission 
towers  was  reported  by  Kawai  [6].  He  used  a  direct 
method  to  measure  tower  surge  impedance.  His  experi¬ 
mental  results  showed  that  the  tower  response  to  a  ver¬ 
tical  current  is  different  from  the  response  to  a  horizon¬ 
tal  current.  Scale-model  measurements  were  reported  by 
Chisholm  [7],  [8]  and  Wahab  et  al.  [9].  These  measure¬ 
ments  results  showed  that  the  tower  surge  impedance  is 
strongly  influenced  by  the  angle  of  current  injection. 

Recently,  theoretical  work  was  reported  by  Ishii  and 
Baba  [10].  They  estimated  the  surge  response  of  a  tower 
by  numerical  electromagnetic  field  analysis.  The  calcu¬ 
lated  results  were  compared  with  the  field  test  results 
[11].  The  analysis  showed  that  surge  response  and  surge 
impedance  of  the  tower  depend  on  the  arrangement  of 


the  current  lead  wire. 

One  of  the  authors  derived  the  formula  of  surge 
impedance  with  ground  surface:  Z  —  60’{ln(2y/2h/ro)- 
1.983}  (11)  and  without  ground  surface:  Z  =  60  • 
{ln(2y/2h/ro)  -  1.540}  (fl)  [12].  The  theoretical  formula 
of  surge  impedance  with  the  ground  surface  is  very  close 
to  the  well  known  experimental  formula  of  Hara  et  al 
[13].  In  this  paper,  we  investigate  the  surge  impedance  of 
the  vertical  conductor  on  the  basis  of  experimental  and 
simulation  analysis.  These  analysis  results  agree  satis¬ 
factorily  with  the  theoretical  values. 

II.  METHOD  OF  ANALYSIS 

For  the  present  analysis,  the  Numerical  Electromag¬ 
netic  Code  (NEC-2)  is  employed.  It  is  a  widely  used 
three-dimensional  electromagnetic  modeling  code  based 
on  the  method  of  moments  [14]  in  the  frequency  domain, 
and  is  particularly  effective  in  analyzing  the  electromag¬ 
netic  response  of  antennas  or  of  other  metallic  structures 
composed  of  thin  wires.  A  vertical  conductor  system 
needs  to  be  decomposed  into  thin  wire  elements,  and  the 
position,  orientation  and  the  radius  of  each  element  con¬ 
stitute  the  input  data,  along  with  the  description  of  the 
source  and  frequencies  to  be  analyzed.  In  the  analysis, 
ail  the  elements  in  the  systems  are  treated  as  perfect 
conductors.  To  solve  the  time-varying  electromagnetic 
response,  Fourier  transform  and  inverse  Fourier  trans¬ 
form  are  used. 

The  validity  of  the  computed  results  when  NEC-2  is 
applied  to  the  analysis  of  surge  response  of  a  vertical  con¬ 
ductor  has  been  verified  by  comparing  with  experimental 
results.  In  the  simulation  and  experimental  analysis,  a 
reduced-scale  model  is  chosen  in  order  to  make  the  ex¬ 
periment  simple  and  flexible.  However,  it  is  not  possible 
to  achieve  the  same  accuracy  as  with  the  full-scale  model, 
especially  in  simulating  the  direct  method,  since  the  ge¬ 
ometrical  size  of  the  measuring  devices  is  large  relative 
to  the  whole  system. 

In  the  lightning  surge  analysis  with  Electromagnetic 
Transients  Program  (EMTP),  the  vertical  conductor 
model  has  been  represented  by  an  equivalent  circuit  of 
the  transmission-line  type  since  it  can  be  easily  inter¬ 
faced  with  EMTP.  Then,  the  basic  parameters  are  its 
surge  impedance,  the  travelling  wave  propagation  veloc¬ 
ity,  and  the  attenuation  and  deformation  characteristics 
of  the  travelling  wave.  Of  these,  the  attenuation  and 
deformation  characteristics  of  the  travelling  wave  deter¬ 
mine  the  reflected  wave  from  the  base  of  vertical  conduc¬ 
tor;  that  is,  the  wave  tail  after  its  peak  in  the  vertical 
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conductor’s  potential-rise  waveform. 

III.  THEORETICAL  FORMULA  OF  SURGE 
IMPEDANCES 

One  of  the  author’s  theory  is  able  to  apply  widely  in 
case  of  ground  surface  and  without  ground  surface.  Sup¬ 
pose  that  the  surge  electric  current  invades  to  the  verti¬ 
cal  conductor  whose  height  is  h  and  radius  is  r0.  Then 
the  surge  current  wave  is  reflected  at  the  ground  of  the 
perfect  conductor  and  returns  to  the  top  of  the  vertical 
conductor. 

Introducing  the  electric  current  reflectivity  0=1  and 
the  magnetic  field  reflectivity  7(7*,  7r)  =  0,  the  theoret¬ 
ical  formula  of  surge  impedance  which  is  very  close  to 
the  well  known  experimental  formula  [13]  is  obtained  as 
follows; 

Z  =  60  •  (ln(^-)  -  b 
v  v2 r0J  A' 


=  60.  (b(^)- 1.983)  (1) 

However,  if  it  is  considered  that  /3  =  7*  =  7r  =  1,  it 
became 

V(t)  =  (,D  (Ct  +  2r°)  _  \ 

U  27T  V  2r0  2  (ct  +  r0)J 

The  above  equation  can  be  modified  by  substituting 
ct  =  2h  and  assuming  h  ro  as  follows; 

z  =  60-(ln(A)_|) 


=  60  .(ln(^)~  1.540)  (2) 

On  the  other  hand,  if  there  is  no  ground,  the  following 
formula  is  induced  [15]. 

V(t)=  f(-ErdI) 

Jo 

__  (ct  +  2rp)  ct  \ 

2n  \n  2 r0  ~  2 (ct  +  r0) ) 

Substituting  ct  =  2h  and  assuming  h  >  r0  in  the 
above  equation,  we  get 

Z  =  60(ln(A)_l) 

To  l 

=  60  (ln(^)- 1.540)  (3) 

This  formula  given  by  (3)  is  the  same  as  (2). 


(a)  With  ground  surface. 


(b)  Without  ground  surface. 


Fig.  1.  Arrangement  of  the  vertical  conductor  system. 


TABLE  I 

Specifications  Of  Measuring  Equipment 


Current 

sensor 

Equipment 

Tektronix 

CT-1 

Frequency 

25kHz~ 

1GHz 

Sensitivity 
5mV/mA  ±  3% 
into  50H  load. 

Voltage 

sensor 

Tektronix 

P6243 

DC~ 

1GHz 

10:1 

Power 

supply 

Tektronix 

1103 

40~ 

440Hz 

±  5VDC 
±  2% 

Recording 

equipment 

HP54616B 

DC~ 

500MHz 

2Gsa/s, 

8  bit  word 

Pulse 

Generator 

HP8131A 

DC~ 

500MHz 

100mV~ 

5Vpp  into  50  ohm 

IV.  EXPERIMENTAL  AND  SIMULATION 
ANALYSIS  OF  SURGE  IMPEDANCE 

Fig.  1  shows  the  reduce-scale  model  of  the  vertical  con¬ 
ductor  system  for  the  simulation  analysis.  The  arrange¬ 
ment  of  the  current  lead  wire  connected  to  the  top  of  the 
vertical  conductor  with  the  existence  of  the  ground  sur¬ 
face  and  without  ground  surface  are  indicated  in  Fig.  1(a) 
and  1(b)  respectively.  Whereas,  Fig.  1(b)  is  also  the  case 
of  the  lightning  phenomena  caused  by  the  return  stroke 
[16].  The  other  type  of  the  lightning  phenomena  caused 
by  a  downward  travelling  current  wave  can  also  be  ex¬ 
amined.  For  the  simulation  of  this  situation,  a  pulse 
current  generator  needs  to  be  placed  remotely  above  the 
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1  ns/div 

(a)  Measured  voltage  waveform. 


1  ns/div 

(b)  Measured  current  waveform. 

Fig.  2.  Experimental  results  of  voltage  and  current  considering 
with  ground  surface. 


channel.  However,  the  computed  waveforms  for  this  case 
of  current  injection,  a  current  wave  travelling  down  the 
current  lead  wire,  are  similar  to  those  of  return  stroke 
type. 

The  experimental  arrangement  of  the  vertical  conduc¬ 
tor  system  with  the  existence  of  ground  surface  is  a  little 
different  from  the  simulation  arrangement  and  is  indi¬ 
cated  in  Fig.  5  of  [15].  However,  the  experimental  ar¬ 
rangement  without  ground  surface  is  similar  to  Fig.  1(b). 
A  voltage  measuring  wire  of  200  cm  in  length  is  placed 
perpendicular  to  the  current  lead  wire  and  is  connected 
to  the  top  of  the  vertical  conductor  which  is  60  cm  in 
height  and  radius  of  0.05  cm.  The  ends  of  the  horizon¬ 
tal  voltage  measuring  wire  in  both  cases  are  stretched 
down  and  connected  to  the  ground  through  matching  re¬ 
sistance.  This  termination  condition  does  not  affect  the 
phenomena  at  the  vertical  conductor  within  17.33  ns. 

A  step  current  pulse  generator  having  pulse  voltage 
of  5  V  in  magnitude,  rise-time  of  1  ns  and  pulse  width 
of  40  ns  is  installed  as  indicated  in  both  cases  which 
is  meant  to  incorporate  the  influence  of  the  induction 
from  the  lightning  channel  hitting  the  vertical  conduc¬ 
tor.  For  the  simulation  analysis,  and  to  save  the  com¬ 
putation  time,  the  conductors  of  the  system  are  divided 
into  10  cm  segments.  To  evaluate  the  voltage  of  the  top 


(a)  Computed  voltage  waveform. 


(b)  Computed  current  waveforms. 


Fig.  3.  Computed  waveforms  of  voltage  at  the  top  and  currents 
in  the  various  parts  of  the  vertical  conductor  in  case  with  ground 
surface. 

of  a  structure,  10  kft  resistance  was  inserted  between  the 
top  of  the  structure  and  the  end  of  the  voltage  measur¬ 
ing  wire.  The  voltage  at  the  top  of  vertical  conductor 
is  measured  by  a  voltage  probe  with  high  resistance  and 
low  capacitance(lMf2  and  jlpF).  The  injection  current 
is  measured  by  current  transformer.  The  specifications 
of  the  measuring  equipment  are  shown  in  Table  I.  The 
waveform  of  current  flowing  through  the  vertical  conduc¬ 
tor  is  also  obtained  from  the  experiment  and  simulation 
analysis.  The  system  of  structures  under  those  analysis 
was  postulated  to  be  on  the  perfectly  conducting  ground. 
Then  we  calculate  the  surge  impedance  which  is  defined 
by  the  ratio  of  the  instantaneous  values  of  the  voltage  to 
the  current  at  the  moment  of  voltage  peak.  * 

As  the  pulse  applied  to  the  current  lead  wire  according 
to  Fig.  1(b),  the  current  starts  flowing  through  the  ver¬ 
tical  conductor  instantly.  However,  for  the  arrangement 
of  Fig.  1(a),  the  current  through  the  vertical  conductor  is 
delayed  by  the  round-trip  time  of  the  travelling  wave  in 
the  conductor.  While  in  both  cases,  the  reflection  wave 
from  the  ground  reaches  the  top  of  the  vertical  conductor 
at  t  —  2h/c,  where  c  is  the  velocity  of  light.  And  that 
is  why  the  maximum  potential  of  the  vertical  conductor 
will  occur  at  time  t  =  2 h/c. 

A.  With  Ground  Surface 

Considering  with  Fig.  1(a),  we  want  to  find  the  volt¬ 
ages  and  currents  experimentally  and  with  the  simula¬ 
tion  by  the  NEC-2.  Fig.  2  shows  the  experimental  re¬ 
sults  of  the  voltage  and  current  with  the  ground  sur- 
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Fig.  4.  Surge  impedances  of  the  vertical  conductor  with  the  ground 
surface  at  0  <  t  <  2 h/c. 


1  ns/div 

(a)  Measured  voltage  waveform. 


1  ns/div 

(b)  Measured  current  waveforms. 


Fig.  5.  Experimental  results  of  voltage  and  current  considering 
with  no  ground  surface. 


face.  Fig.  3  shows  the  simulation  results  by  the  NEC- 
2.  In  Fig.  3(a),  the  influence  of  the  reflected  wave  from 
the  ground  reaches  the  top  of  the  conductor  is  observed 
at  t  —  2 h/c  =  4  ns  exactly  which  means  that  the 
travelling  wave  is  propagating  at  the  velocity  of  light. 
Fig.  3(b)  shows  the  computed  waveforms  of  current  flow¬ 
ing  through  the  vertical  conductor  as  indicated  by  the 
mark  ‘CT?  in  Fig.  1(a).  As  the  pulse  generator  is  placed 
300  cm  from  the  vertical  conductor,  the  current  through 
the  vertical  conductor  is  being  delayed  approximately  10 
ns.  The  waveforms  start  rising  after  10  ns  which  can  be 
noticed  from  Fig.  3.  The  existence  of  the  ground  surface 
can  be  observed  in  Fig.  3(b).  where  the  field  produced 


Fig.  6.  Computed  waveforms  of  voltage  at  the  top  and  currents  in 
the  various  parts  of  the  vertical  conductor  in  case  with  no  ground 
surface. 


Fig.  7.  Surge  impedances  of  the  vertical  conductor  at  0  <  t  <  2h/c 
with  no  ground  surface. 


by  the  current  injected  horizontally  induce  current  of 
small  magnitude  before  the  actual  surge  current  flowing 
through  the  vertical  conductor.  These  simulation  results 
of  currents  in  Fig.  3(b)  obtained  by  the  NEC-2  exactly 
coincide  with  the  experimental  results  [15].  Then,  we 
compare  the  theoretical  value  of  the  surge  impedance 
considering  the  ground  surface  given  by  the  (1)  with  the 
simulation  and  experimental  results  of  that.  Fig.  4  shows 
that  the  vertical  conductor  surge  impedances.  The  theo¬ 
retical  values  of  surge  impedance  calculated  by  using  (1) 
is  just  after  the  surge  electric  current  reaches  the  ground 
and  produce  reflected  current  wave.  As  we  need  to  know 
surge  impedance  at  t  =  2 h/c  —  4  ns.  In  these  results. 
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(a)  With  ground  surface. 


5V,  Ins 

pulse  # 
generator 


J  500 

(b)  Without  ground  surface. 


Fig.  8.  Equivalent  circuits  of  the  vertical  conductor  models. 


theoretical,  computed  and  experimental  values  of  surge 
impedances  are  approaching  closely  at  2h/c. 

B.  Without  Ground  Surface 

Fig.  5  shows  the  experimental  results  of  the  voltage 
across  the  voltage  measuring  wire  and  currents  through 
upper  and  lower  parts  of  the  vertical  conductor  in  the  ab¬ 
sence  of  ground  surface.  The  simulation  results  of  volt¬ 
age  at  the  top  of  the  vertical  conductor  and  currents 
through  different  parts  of  it  are  shown  in  Fig.  6.  How¬ 
ever,  in  this  case  of  analysis,  the  waveforms  of  current 
through  the  vertical  conductor  are  somewhat  different 
from  Fig.  3(b)  at  the  starting  region  because  of  absence 
of  the  ground  surface.  Also,  the  current  starts  flowing 
instantly  through  the  vertical  conductor  without  being 
delayed.  Finally,  the  theoretical,  computed  and  mea¬ 
sured  values  of  surge  impedances  are  shown  in  Fig.  7. 
Here  also  we  see  that  the  values  of  surge  impedances 
approach  closely  at  t  »  2 h/c. 

V.  VERTICAL  CONDUCTOR  MODELS  FOR 
EMTP  ANALYSIS 

The  Electromagnetic  Transients  Program  (EMTP)  is 
probably  the  most  widely-used  power  system  transients 


(a)  Voltages  across  the  voltage  measuring  resistance. 


(b)  Currents  through  the  vertical  conductor. 


Fig.  9.  Simulation  results  of  voltages  and  currents  by  the  EMTP 
at  0  <  t  <  2 h/c. 


Fig.  10.  Simulation  results  of  surge  impedances  by  the  EMTP  at 
0  <  t  <  2 h/c. 


simulation  programs  in  the  world  today.  In  this  sec¬ 
tion,  the  EMTP  simulations  based  on  the  circuit  the¬ 
ory  were  perfomed  for  the  vertical  conductor  model  with 
ground  and  without  ground  surface.  In  the  circuit  model, 
the  line  was  represented  by  a  distributed  R-L-C  circuit 
with  the  skin  effect  being  neglected.  The  NEC-2  cannot 
exactly  model  the  structures  of  actual  towers  or  verti¬ 
cal  conductors,  and  in  addition  it  cannot  directly  inter¬ 
faced  with  the  EMTP.  For  the  EMTP  simulations,  there¬ 
fore,  it  is  practical  to  employ  an  equivalent  circuit  of  the 
transmission-line  type  for  representing  the  vertical  con¬ 
ductor  system.  In  developing  the  model  or  in  determin¬ 
ing  its  parameters,  characteristics  stated  in  the  preceding 
sections  should  be  taken  into  consideration.  In  this  sec¬ 
tion,  vertical  conductor  models  used  so  far  are  reviewed 
with  emphasis  on  their  performance  in  reproduction  of 
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measured  waveforms  of  current  through  the  vertical  con¬ 
ductor  and  voltage  at  the  top  of  it.  The  surge  impedance 
is  then  calculated  from  the  ratio  of  the  maximum  poten¬ 
tial  at  the  conductor  top  to  the  current  through  it  at  the 
time  of  voltage  peak.  Fig.  8  shows  the  equivalent  circuit 
representation  for  the  vertical  conductor  system  consid¬ 
ering  with  the  existence  of  ground  surface  and  without 
ground  surface.  The  dimensions  of  these  circuit  models 
are  the  same  as  considered  for  the  simulation  and  experi¬ 
mental  model  systems  of  Fig.  1.  The  voltage  sensors  and 
the  current  sensors  indicated  in  Fig.  8  represent  the  mea¬ 
suring  points.  The  surge  impedances  of  the  distributed 
line  are  used  as  the  input  data  for  the  EMTP  analysis, 
and  are  of  different  values  depending  on  the  height  of 
the  conductors.  As  the  analysis  with  EMTP,  it  can  be 
easily  handled  to  the  horizontal  conductor  but  cannot 
be  handled  just  as  it  is  to  the  perpendicular  conductor. 
Therefore,  to  solve  the  problem,  the  perpendicular  con¬ 
ductor  can  be  divided  into  the  horizontal  conductors  as 
it  makes  a  center  level  at  the  axis  in  each  conductor. 

Fig.  9  shows  the  simulation  results  of  voltages  and  cur¬ 
rents  with  the  existence  of  ground  surface  and  without 
the  ground  surface  by  the  EMTP  analysis  for  the  equiv¬ 
alent  circuit  representation  of  Fig.  8.  The  solid  lines  of 
Fig.  9  correspond  to  circuit  representation  of  Fig.  8(b) 
and  the  chain  lines  for  the  Fig.  8(a).  The  starting  time  of 
current  flowing  through  the  vertical  conductor  depends 
on  the  position  of  the  pulse  generator.  As  the  pulse  is  in¬ 
jected  at  300  cm  from  the  vertical  conductor  with  ground 
surface  as  in  Fig.  8(a),  the  currents  start  flowing  after  10 
ns  of  currents  through  it  in  case  with  no  ground  surface. 
The  occurrence  of  the  reflection  can  also  be  observed  in 
Fig.  9.  Then  the  EMTP  results  of  the  surge  impedances 
of  the  vertical  conductor  with  ground  surface  and  with¬ 
out  ground  surface  are  shown  in  Fig.  10. 

TABLE  II 

Surge  Impedances  of  the  Vertical  Conductor  at  t  ss  2h/c 


With  ground 

Without  ground 

Theoretical 

368 

395 

Computed 

373 

402 

Measured 

366 

406 

EMTP 

375 

403 

Theoretical,  experimental  and  computed  by  the  NEC- 
2  and  EMTP  results  can  be  summarized  in  the  Table 
II  at  t  =  2h/c  so  as  to  make  quantitative  evaluation. 
The  surge  impedance  for  the  ground  surface  is  naturally 
much  lower  than  without  ground  surface  that  can  also 
be  realized  by  the  (1)  and  (2).  The  theoretical  values 
of  surge  impedance  agree  well  with  the  computed  and 
experimental  values. 


VI.  CONCLUSIONS 

The  theoretical  values  of  surge  impedances  are  verified 
by  comparing  the  computed  and  experimental  results  on 
simple  structures.  The  difference  is  less  than  about  5%, 
which  is  within  the  accuracy  maintained  in  the  analysis. 
Also,  the  travelling  wave  propagates  at  nearly  the  veloc¬ 
ity  of  light.  The  surge  characterist  ics  have  some  influence 
on  the  type  of  the  lightning  current  with  the  presence  of 
ground  surface  and  without  the  ground  surface.  The 
difference  comes  from  the  different  electromagnetic  field 
around  the  vertical  conductor  influenced  mainly  by  the 
electric  fields  associated  with  the  currents  propagating 
the  vertical  conductor  and  current  lead  wire.  Also  the 
restriction  of  the  size  of  the  perfectly  conducting  ground 
plane  and  the  effect  of  the  voltage  probes  might  cause 
small  difference  in  the  experimental  results  of  voltage 
and  current  waveforms. 
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ABSTRACT 

At  present  multi-GHz  operating  frequencies,  the  elec¬ 
trical  properties  of  high-end,  multilevel  IC  intercon¬ 
nects  must  be  described  with  Maxwell’s  equations. 
We  have  developed  an  entirely  new  floating  random- 
walk  (RW)  algorithm  to  solve  the  2D  time-harmonic 
Maxwell-Helmholtz  equation.  The  algorithm  requires 
no  numerical  mesh,  thus  consuming  a  minimum  of 
computational  memory — even  in  complicated  prob¬ 
lem  domains,  such  as  those  encountered  in  IC  inter¬ 
connects.  The  major  theoretical  challenge  of  deriving 
an  analytical  Green’s  functions  in  arbitrary  heteroge¬ 
neous  problem  domains  has  been  successfully  re¬ 
solved  by  means  of  an  accurate  approximation:  itera¬ 
tive  perturbation  theory.  Initial  numerical  verification 
of  the  algorithm  has  been  achieved  for  the  case  of  a 
“skin-effect”  problem  within  a  uniform  circular  con¬ 
ductor  cross  section,  and  also  for  a  heterogeneous 
“split-conductor”  problem,  where  one  segment  of  a 
square  domain  is  conducting  material,  while  the  other 
segment  is  insulating.  As  an  example  of  electrical 
parameter  extraction  using  this  algorithm,  we  have 
extracted  the  frequency-dependent  impedance  of  the 
uniform  circular  cross-section  previously  mentioned. 
Excellent  agreement  has  been  obtained  between  the 
analytical  and  RW  solutions,  supporting  the  theoreti¬ 
cal  formulation  presented  here. 

Index  Terms — Floating  random-walk,  Helmholtz 
equation,  Maxwell  equations,  perturbation  theory, 
skin  effect,  IC  interconnect. 


L  INTRODUCTION 

Advances  in  digital  IC  technology  have  resulted  in 
multi-GHz  operation  frequencies.  At  such  frequen¬ 
cies,  circuit  designers  must  account  for  electromag¬ 
netic  phenomena  that  are  difficult  to  calculate.  They 
include  skin-effect  loss,  frequency-dependent  induc¬ 
tance  and  capacitance,  slow-wave  substrate  coupling, 
distributed  transmission-line  propagation  and  high- 
frequency  radiation.  Our  principal  objective  here  is  to 
invent  a  new  numerical  algorithm  capable  of  effi¬ 
ciently  describing  these  increasingly  significant  elec¬ 
tromagnetic  phenomena.  Our  hope  is  to  establish  a 
new  approach  for  the  modeling  and  design  of  com¬ 
plex,  multilevel  IC-interconnects. 

Traditional  numerical  methods  for  solving  electro¬ 
magnetic  problems,  unfortunately,  require  a  discreti¬ 
zation  mesh.  Mesh  size  and  the  resultant  difficulty  of 
solution  become  somewhat  unmanageable  in  compli¬ 
cated  3D  problem  domains.  The  random  walk  (RW) 
algorithm  [  1  ]  that  we  present  here  does  not  employ  a 
mesh.  In  essence,  the  algorithm  executes  a  Monte 
Carlo  integration  [2]  of  an  infinite  series  of  multi¬ 
dimensional  integrals  by  means  of  RWs  through  the 
problem  domain.  These  integrals  contain  “surface” 
and  “volume”  Green’s  function  kernels.  Note,  impor¬ 
tantly,  the  RW  method  is  inherently  parallel,  requir¬ 
ing  minimal  inter-processor  communication. 

A  large  portion  of  the  traditional  RW  literature  treats 
the  Helmholtz  equation  in  homogeneous  problem 
domains  [3].  This  is  principally  because  of  the  ab- 
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sence  of  an  exact  analytical  Green’s  function  in  arbi¬ 
trary  heterogeneous  domains  [4].  The  RW  algorithm 
we  present  in  this  work,  on  the  other  hand  is  applica¬ 
ble  to  heterogeneous  problem  domains — essential  for 
IC-interconnect  modeling. 

The  primary  objective  of  this  work  is  the  detailed 
theoretical  formulation  of  a  novel  floating  RW  algo¬ 
rithm  based  on  iterative  perturbation  theory.  In  Sec¬ 
tion  II,  we  develop  a  vector-potential  formulation  of 
the  2D  Maxwell-Helmholtz  equation,  suitable  for 
skin-effect  analysis.  A  derivation  of  the  relevant 
Green’s  functions  for  the  2D  Maxwell-Helmholtz 
equation  using  iterative  perturbation  theory  is  given 
in  Section  III.  In  Section  IV,  we  apply  the  Green’s 
functions  defined  in  Section  HI  to  define  a  specific 
floating  RW  algorithm.  Section  V  presents  the  results 
of  a  numerical  2D  skin-effect  problem  analysis 
within  a  circular  conductor  cross  section,  including 
the  frequency-dependence  impedance  per  unit  length 
for  the  circular  cross  section  at  different  frequencies. 
This  section  also  contains  the  results  for  a  heteroge¬ 
neous  “split-conductor”  problem,  where  one  segment 
of  a  square  cross  section  is  electrically  conducting, 
while  the  other  is  insulating.  For  each  one  of  these 
problems,  comparison  with  an  exact,  analytical  solu¬ 
tion  is  provided.  Lastly,  Section  VI  summarizes  our 
work  and  indicates  possible  future  directions. 

II.  PROBLEM  FORMULATION 

Consider  a  2D  solid-conductor  cross  section  in  the  xy 
plane,  where  we  impress  a  z-directed  current  density 
at  the  conductor  surface.  We  define  a  corresponding 
current-density  phasor  Jz  in  the  harmonic  steady  state. 
We,  furthermore,  neglect  any  free-charge  density  as 
an  approximation.  Time-harmonic  Maxwell’s  equa¬ 
tions  require  that  the  electric-field  phasor  within  the 
conductor  cross  section  satisfy  the  scalar  Helmholtz 
equation  [5]: 

V2Ez-y2E.=  0.  (1) 

Above,  V2  =  d2  Idx2  +d2  ldy2\  E.  =E.(x,y)-, 
y2  =  -£i0£6)2  +ijjQ<7(o  ;  jUo,  e  and  crare  operation 

frequency,  free-space  magnetic  permeability,  permit¬ 
tivity  and  conductivity,  respectively.  At  the  conduc¬ 
tor  surface,  the  impressed  current  density  expresses 
itself  as  a  boundary  condition  in  electric  field  by 
means  of  the  Ohm’s  Law  constitutive  relation: 


£z=  — •  (2) 

a 

Equations  (1)  and  (2)  essentially  describe  the  so- 
called  2D  “skin-effect  problem”  in  our  conductor. 
Electric  field,  or  equivalently,  current  density,  will 
vary  within  the  conductor  as  a  function  of  frequency 
and  material  parameters,  subject  to  an  applied  surface 
boundary  condition.  We  choose  now  to  reformulate 
the  problem,  using  vector  potential  A  =  Az(x,  y)e., 
with  V  x  A  =  B  in  the  Coulomb  gauge  V  •  A  =  0  [6]. 
This  formulation  is  useful  in  a  future  3D  extension  of 
this  work,  because  it  conveniently  decouples  field 
components  in  the  governing  equations. 

Equations  (1)  and  (2),  in  the  vector-potential  formu¬ 
lation,  generate  a  “forced”  Maxwell-Helmholtz  sys¬ 
tem: 

V2Az-y2Az=-/i0<r^,  (3) 

di 

Jz  =-icroA.  -<7-^,  (4) 


where,  at  the  conductor  surface, 

A.  =0.  (5) 

The  quantity  - icroA z  above  is  the  so-called  “eddy- 
current  density”.  In  deriving  (3)  and  (4)  from  Max¬ 
well’s  equations,  we  observe  that  for  no  free-charge 
density,  the  scalar  potential  function  <p  is  frequency 
independent  and  it  is  completely  decoupled  from 
vector  potential  Az. 

In  addition,  as  d(pldi  generally  depends  solely  on  x 
and  y,  and  not  ox  it  must,  as  well,  satisfy  (4)  in  the  dc 
limit  <y-4  0.  Accordingly,  -ad(pl  di  can  be  identi¬ 
fied  as  the  dc  current  density  phasor.  It  should  be 
noted,  that  though  this  phasor  has  a  non-constant 
harmonic  temporal  variation  exp(/&r)  for  any  a>±  0, 
its  spatial  dependence  remains  identical  to  that  at  dc. 

We  require,  as  well,  correspondence  with  surface 
condition  (2).  We  must  be  careful,  therefore,  to  im¬ 
press  a  surface  current  density  in  (2)  consistent  with 
(4)  and  (5).  In  other  words,  we  define  an  ideal  source 
as  one  that  excites  our  conductor  cross  section  ac¬ 
cording  to  (4)  and  (5).  An  ideal  source  maintains  the 
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dc  current-density  phasor  spatial  dependence  at  the 
conductor  surface,  with  the  proviso  that  the  phasor 
temporal  dependence  is  of  the  form  exp(ier). 

ffl.  ITERATIVE  PERTURBATION  THORY 
BASED  GREEN’S  FUNCTION 

The  Green’s  function  equation  corresponding  to  the 
2D  Maxwell-Helmholtz  equation  (3)  is 

W2G-y2G  =  S(  r-r0),  (6) 

where  G  =  G(r  I  r0)  is  the  Green’s  function  at  (x,  y ) 
position  coordinate  r  due  to  a  2D  Dirac  delta- 
function  source  at  r0.  Equation  (6)  does  not,  gener¬ 
ally,  have  an  analytical  solution  for  arbitrary  }(r). 
We  will  derive,  nonetheless,  using  iterative  perturba¬ 
tion  theory  [7],  an  approximate  expression  for  G  on 
the  circular  domain,  with  arbitrary  radius  R,  shown  in 
Fig.  1.  This  Green’s  function  will  allow  us  to  later 
develop  a  novel  RW  algorithm  for  the  solving  2D 
skin-effect  equation  (3).  The  Green’s  function  G  is 
assumed  to  be  zero  on  the  boundary  of  the  circular 
domain,  as  the  problems  under  consideration  are 
Dirichlet  [8]  problems. 


Figure  1:  A  circle  of  arbitrary  radius  R  over  which  the 
Green’s  function  in  (6)  is  estimated. 


Let  us  define  the  zeroth-order  approximation  G(m  for 
G  as  the  solution  to  (6)  with  y=  0.  Therefore, 

V2G(0)  =£(r-r0).  (7) 

Above,  r  (p,  6)  is  the  point  where  the  zeroth-order 
approximation  is  calculated  given  a  delta  function 
centered  at  r0  (p0,  d0). 


Using  (6)  for  iteration,  we  can  then  generate  a  first- 
order  approximation  CjU  in  terms  of  G0): 

V2G(I)  =S(r-ro)  +  y2G(0).  (8) 

The  solution  to  Poisson  equation  (7)  is  well  known;  it 
has  the  form,  in  polar  coordinates  [9] 

G(0)  =  — lnl"/?2  — 

4 7t  L  B 

A  =  p2  +pl-2pp0cos(e-eo).  (9) 

B  =  p2p20  +R4 -2pp0R2  cos(0-0o). 

Now,  we  are  in  a  position  to  evaluate  G(,)  from  (8). 
Using  the  expression  for  G(0),  and  with  the  right  side 
of  (8)  as  the  Poisson  source  term,  we  find  an  expres¬ 
sion  for  the  first-order  approximation  to  (6)  given  by 
[10] 

G(,)(rlr0)  =  JJ</2r5G(0)  (rlrs) 

[<J(rs  - r0 )  +  y2 (rs  )G(0)  (r5  I  r0 )  (1Q) 

=  G(0)  (r  lr0)  + 

P%X2(rs)G(0)(rlrs)G(0)(rslro). 

s 

Note  that  G(,)  given  by  (10),  is  an  approximate  ex¬ 
pression  for  G  as  given  by  (6).  The  integration  vari¬ 
able  in  (10)  represents  an  infinitesimal  area  element 
on  the  circular-domain  surface  S  in  Fig.  1.  Note,  as 
mentioned  earlier,  homogeneous  Dirichlet  conditions 
have  been  employed  in  obtaining  (9)  and  (10). 

We  next  use  this  approximate  Green’s  function  G^'to 
develop  a  general  solution  to  skin-effect  equation  (3) 
within  our  circular  domain  in  Fig.  1.  Two  integral 
terms  arise — a  line  integral  about  the  domain  circum¬ 
ference  C  which  takes  into  account  the  effect  of 
boundary  conditions,  and  a  surface  integral  through¬ 
out  the  domain  S  itself,  which  takes  into  account  the 
effect  of  the  source  term,  and  the  vector-potential  at 
the  center  of  the  circular  domain  is  given  by  [1 1] 
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Az  (center)  =  <fAz  (/?,  &)VTq  G(r  I  r0  >n  dc 

c 


-Mo<* 


(rlr 0)ds. 


(11) 


Substituting  (10)  in  (1 1)  and  after  some  mathematical 
manipulation,  we  obtain,  for  A,  at  the  domain  center 


A.  (center) «  ft* dd  Az (R,  6) 


1  |  1 
2  7C  4 n1 


i  wqm 
q 


/ 

V 


\ 

J 


ftdpft*ddp 


1 

8  K1 


I  Fq(p,0) 
<? 


(12) 


where 


C  =  Tjln(rj/R)(R2-tj2)  (13) 

D  =  R2  +T)2  -2Rrjcos(0-^). 


and 


Fq (A *)  =  Y]  ft dT}^_,2K/2d£ tj ln(tj /  R) In  — 

E  =  p2R2  +jj2R2  -2prjR2  cos(0-g)  (14) 
F  =  p2 Tj2  +R*  -2prjR2  cos(0-g). 


For  simplicity,  above,  we  take  f  to  be  piecewise  con- 
stant  with  respective  values  yq  in  0 -quadrants  q  =  1, 
2,  3,  and  4.  The  quantities  within  square  brackets  in 
(12)  are  2D  versions,  respectively,  of  surface  and 
volume  Green’s  functions  encountered  in  3D  prob¬ 
lem  domains.  These  Green’s  functions  consist  of  two 
auxiliary  functions  Wq  and  Fq  defined  in  (13)  and 
(14).  The  functions  represent  perturbative  corrections 
arising  from  the  fAz  term  in  the  original  Maxwell- 
Helmholtz  equation  (3).  In  (13)  and  (14),  ^and  £are 
variables  of  integration.  7]  takes  values  between  0 


and  /?,  while  6  assumes  values  between  {q-\)nf2 
to  qn  1 2  for  a  particular  quadrant.  Equations  (12)- 

(14)  are  the  starting  point  for  defining  a  RW  algo¬ 
rithm  for  solving  (3)  in  2D  domains  with  arbitrary 
piecewise-constant  spatial  variation  in  y  subject  to 
arbitrary  Dirichlet  boundary  conditions. 

The  total  current,  /,  through  the  cross  section  can  be 
calculated  by  integrating  the  current  density  given  in 
(4)  over  the  problem  domain  ( ds  being  an  infinitesi¬ 
mal  area  unit)  and  can  be  written  as 


/  =  ]\ds 
s 


-iaccAz 


(15) 


The  integral  expression  for  vector  potential  from  (12) 
is  substituted  in  (15)  to  obtain  a  multi-dimensional 
integral  expression  for  total  current  through  the  con¬ 
ductor  surface. 

The  internal  impedance  per  unit  length  is  defined 
as[12] 


E  .(devalue)  =  — 

Z,=— - - - —  •  06) 

At  this  point,  the  crucial  thing  to  note  is  that  for  esti¬ 
mating  frequency-dependent  impedance,  we  need  not 
estimate  field  or  vector  potential  at  a  large  number  of 
points  within  the  problem  domain,  or  for  that  matter, 
at  any  point  within  the  problem  domain.  The  problem 
of  impedance  extraction  is  reduced  to  estimating  the 
overall  multi-dimensional  integral  expression  for 
current  obtained  from  (15)  within  the  FRW  frame¬ 
work  to  be  described  in  the  next  section,  and  then 
using  (16)  to  evaluate  the  internal  impedance  per  unit 
length. 


IV.  THE  FLOATING  RW  ALGORITHM 

As  mentioned  earlier,  the  floating  RW  algorithm  is  a 
Monte  Carlo  evaluation  of  an  infinite  series  of  multi¬ 
dimensional  integrals.  The  kernels  of  these  integrals 
consist  of  products  of  surface  and  volume  Green’s 
functions.  In  this  section,  we  describe  the  floating 
RW  algorithm  in  detail  in  context  of  the  skin-effect 
problem  in  a  circular  cross  section.  As  shown  in  Fig. 
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2,  we  define  RWs  to  start  at  a  point,  where  we  need  ^0)  _  ^  ^  ^ 

to  estimate  Az  in  (3). 


Figure  2:  A  schematic  diagram  of  a  circular  cross  section  is 
shown.  One-,  two-  and  a  three-  hop  RWs  are  represented. 

The  RWs  propagate  as  “hops”  of  different  sizes  from 
circle  centers  to  circumferences,  consistent  with  a 
statistical  interpretation  [1]  of  (12).  Maximally  sized 
circles,  subject  to  limitations  imposed  by  iterative 
perturbation  theory,  are  used  with  hop-location  prob¬ 
ability  rules  again  consistent  with  (12). 

We  define,  with  each  hop,  a  numerical  weight  factor 
derived  from  (12).  The  product  of  these  weight  fac¬ 
tors  over  a  walk,  multiplied  by  the  solution  at  the 
problem  boundary — where  the  walk  must  termi¬ 
nate — gives  a  statistical  estimate  for  Az  at  the  RW 

starting  point.  All  this,  again,  is  entirely  consistent 
with  a  statistical  interpretation  of  (12).  We  can  thus 
obtain  an  accurate  statistical  estimate  for  Az  by  av¬ 
eraging  over  a  statistically  large  number  of  RWs. 
Mathematically,  we  can  write  such  an  estimate 

(17) 

N  n= 1 


where  N  is  the  number  of  walks  and  A  ‘M)  is  the  nth- 

walk  estimate.  Referring,  again,  to  Fig.  (2),  we  see 
examples  of  three  representative  RWs:  a  one-hop,  a 
two-hop,  and  a  three-hop  walk.  The  contributions 
from  these  three  RWs  can  be  written  as 


A™  =Ks(rl)+Kc(rl)Ks(r2), 

A^  =  Ks(ri)+Kc(rl)Ks(r2)  + 
Kc  (r,  )KC  (r2  )KS  (r3 ). 


Above,  Kc  represents  the  weight  factor  associated 
with  the  “surface”  Green’s  function,  the  6- 
integral  term  in  (12).  The  function  Ks  represents  the 
weight  factor  associated  with  the  “volume”  Green’s 
function,  the  (/?,$- integral  term  in  (11).  Assuming 
the  hops  are  uniformly  distributed  in  (/?,0),  these 
weight  factors  have  the  form,  from  (12), 


Kc 


=  R 


1 

1  +  — X  WJ0) 
2lt  q  q  J 


(19) 


and 


Ks 


±ln(p/R)+±ZFq(p,0) 

JL  o7T  q 


(20) 


For  estimating  the  frequency-dependent  impedance 
per  unit  length  as  given  in  (16)  a  similar  exercise  is 
carried  out  using  a  statistical  interpretation  of  (15). 
For  heterogeneous  problems,  there  are  a  couple  of 
differences.  First  of  all,  the  maximum  hop  size, 
which  is  decided  by  the  validity  of  iterative  perturba¬ 
tion  theory,  is  different  for  different  medium.  In  this 
paper,  the  maximum  hop  size  is  estimated  to  be  the 
minimum  of  two  numbers.  First,  we  allow  the  first- 
order  correction  in  the  expression  for  the  volumetric 
Green’s  function  given  in  (12)  to  be  equal  to  ten  per¬ 
cent  of  the  zeroth-order  approximation  and  calculate 
a  maximum  hop  size  under  this  assumption.  A  similar 
process  is  carried  out  for  the  surface  Green’s  function 
term  in  (12)  and  a  maximum  hop  size  is  calculated 
under  this  assumption.  The  maximum  hop  size  for 
our  RW  algorithm  is  the  smaller  of  these  two  num¬ 
bers.  Secondly,  the  random  hops  are  restricted  by 
material  interfaces  in  heterogeneous  problems. 
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We  close  this  section  with  a  pseudo-code  listing  that 
defines  our  floating  RW  algorithm  for  estimating 
vector-potential. 

— Floating  Random-Walk  Algorithm  Pseudo- 
Code — 


1)  Choose  the  point  where  Az  need  be  esti¬ 
mated;  call  it  AZ. 

i 

2)  Evaluate  8  -  the  maximum  hop  size  as  de¬ 
termined  by  validity  of  perturbation  theory, 
according  to  the  procedure  described  previ¬ 
ously  in  this  section. 

3)  A- a  pre-defined  small  number. 

4)  NMAX  =  a  pre-defined  large  integer. 

5)  N  =  0. 

6)  TOTALSUM  =  0. 

7)  SUM  =  0. 

8)  Evaluate  the  maximal  radius  that  contacts 
the  closest  problem-domain  boundary ,  with¬ 
out  passing  through  it;  call  it  RMAX. 

9)  RAD  =  MIN  (RMAX,  S  ). 

10)  Draw  a  circle  of  radius  RAD. 

11)  Hop  to  a  point  on  the  circumference  in  con¬ 
formity  with  a  uniform  probability  distribu¬ 
tion  in  6. 

12)  Evaluate  the  exact  weight  factor  Kcfrom 

(13)  and  (19);  call  it  KC. 

13)  Evaluate  the  exact  weight  factor  Ksfrom 

(14)  and  (20);  call  it  KS. 

14)  KC  (zeroth  hop)  =  1;  INCREMENT  =  KC 
( previous  hop)  *  KS  (present  hop). 

15)  SUM  =  SUM  +  INCREMENT. 

16)  IF  (a  boundary  is  not  reached)  THEN  (re¬ 
peat  steps  8  -75). 

17)  IF  (a  boundary  is  reached)  THEN  (termi¬ 
nate  walk ;  N  =  N  +  1 ;  SUM.TOT AL  = 
SUM_TOTAL  +  SUM). 

18)  IF  (N  <  NMAX)  THEN  (repeat  Steps  7-17). 

19)  IF  (N  >  NMAX)  THEN  (AZ  = 
SUM_TOTAL  /  NMAX). 

20)  Evaluate  exact ,  analytical  solu¬ 
tion  AZ(exact) . 

21)  ERROR  =  |AZ  -  AZ(exact)| . 

22)  IF  (ERROR  >  A)  THEN  (NMAX  =  NMAX 
*  1.2;  repeat  steps  5  -21). 

23)  IF  (ERROR  <  A)  THEN  (AZ  =  estimated 

value  ofAr). _ 


V.  VERIFICATION  WITH  THE  HELP  OF 
BENCHMARK  PROBLEMS 

The  principal  objective  of  this  work  is  to  formulate 
and  to  define  a  novel  RW  algorithm  for  2D  Maxwell- 
Helmholtz  equation  solution.  We  have  benchmarked 
our  formulation  against  two  known  solutions.  As  said 
earlier,  the  first  problem  is  a  single  circular  cross 
section,  where  an  alternating  current  of  single  fre¬ 
quency  is  impressed.  Using  the  algorithm  developed 
earlier,  we  estimate  the  current  density  profile  across 
the  cross  section  as  well  as  the  internal  impedance. 
The  analytical  solution  for  the  current  density  along  a 
uniform,  circular-conductor  cross  section  of  radius  R 
is  [12] 


(21) 

Jz(  0) 

where  J0  is  the  zeroth-order  Bessel  function.  The 

variable  p  here  denotes  radial  coordinate  from  the 
conductor  center.  For  this  circular  cross  section,  an 
analytical  expression  for  the  internal  impedance  per 
unit  length  as  defined  in  (16)  is  given  by  [12] 


iyJ0(i}R) 

2nRcrJ0(i}R) 


(22) 


As  mentioned  earlier,  the  second  problem  solved  is  a 
heterogeneous  “split-conductor”  problem,  where  a 
square  domain  is  divided  into  two  unequal  rectangu¬ 
lar  domains  of  insulating  and  conducting  material,  as 
shown  in  Fig.  (3). 

f  -l GHz 


„  l - I 

2  pm\ - 1 

fe  €r  =  2.7 

m  p- 1.8  pCl  •  cm 


Figure  3:  2D,  split-conductor  problem,  the  geometry. 
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The  insulating  region  is  represented  in  light  gray, 
while  the  conducting  region  is  represented  in  dark 
gray.  The  boundary  regions,  where  the  z-component 
of  the  vector  potential  is  known,  are  represented  in 
black.  The  boundary  conditions  are  chosen  such  that 
Az  =  0  in  the  top,  bottom  and  right  boundary  re¬ 
gions,  while  A.  =  sin(/ry  /  L)  in  the  left  boundary 
region  (L  being  the  length  of  the  side  of  the  square) 
where  the  origin  coincides  with  the  left  and  bottom 
comer.  Assuming  the  continuity  of  the  solution  and 
its  derivative  at  the  material  interface  (Lo  being  the 
length  of  the  dielectric),  the  analytical  solution  to  this 
heterogeneous  problem  is  given  by  [13] 


A,  =  sin(Jry) 

[A  sinh(krf  *)+  cosh(fcrf  *)],  0  <  x  <  L0 . 
A,  =sin(Ay) 

[C  sinh{Arc  u  -  d}}  L0  <  x  <  L. 


(23) 


The  constants  in  (23)  are  given  in  (24)  and  (25)  as: 

*  =  *. 

L 

kd  =V^+*2- 

kc=Jrf+k*. 


ay 

o2 

h 

b2 

fli 

a2 

bi 

b2 

Oi 

a3 

by 

0| 

Cl  2 

by 

(24) 


a,  =sinh  (kdL0\a2  =  -sinh[fcf  (L0 -L)l 
a}  =-cosh  {kdL). 

fc,  =ifea],b2  = 

Yd  =-fJ£0)2,Yc 

We  coded  the  algorithm  in  MATLAB  5.0™,  using  a 
400-MHz  Apple  PowerBook  G3™  development  plat¬ 
form.  The  resistivity  for  conducting  material  is  given 
by  p  =  1.8//Q-cm  and  relative  permittivity  of  di- 


.  — 1 — a-.  A  =  Vfflv 
V  to 


(25) 


electric  material  is  given  by  er  =  2.7.  The  radius  of 

the  circular  cross  section  is  given  by  R  -  5//m  while 
for  the  split-conductor  problem,  the  dimensions  of 
the  square  cross  section  is  given  by 
10//mxl0/zm.The  respective  dimensions  of  the  in¬ 
sulating  and  the  conducting  materials  in  the  split- 
conductor  problems  are  shown  in  Fig.  (3).  The  oper¬ 
ating  frequency  /  =  dlK  =  1GHz  corresponds  to  a 
skin  depth  £  =  2.1//m  and  a  wavelength  of 

1.8xl05/fln.The  propagation-constant  squared  (/) 
is  equal  to  4.386x1 0n//m2  within  the  conductor  and 
equal  to  —1185.431197m'2  within  the  dielectric  at  1 
GHz.  Based  on  these  numbers  and  the  criterion  given 
in  Section  IV,  the  maximum  radius  of  hops  inside 
conducting  material  is  0.95//m,  and  the  maximum 

radii  of  hops  within  the  dielectric  material  is 
4 

1.8x10  //m,  which  is  about  twice  the  dimensions  of 

a  chip  (based  on  a  1  cmxlcm  chip).  Thus  we  see  that 
this  perturbation  theory  based  approach  has  the  po¬ 
tential  to  allow  meaningful  interconnect  analysis. 


Figures  (4)  and  (5)  show  the  magnitude  ratio  and 
phase  lag,  respectively,  of  the  skin-effect  current  den¬ 
sity  phasor.  A  total  of  20,000  RWs  were  performed 
per  solution  point.  The  figures  show  excellent  agree¬ 
ment  between  the  analytical  and  RW  solutions.  The 
mean  absolute  error  between  exact  and  RW  solutions 
was  0.001  for  magnitude  and  0.012rad  for  phase. 
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Figure  4:  2D  skin-effect  problem  in  a  homogeneous 
circular  conductor  cross  section,  relative  magnitude. 
Problem  radius  R  =  5//m  and  (o=  1GHz. 
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arg[/»/7;(0)]  (rad) 


p(m) 


Figure  5:  2D  skin-effect  problem  in  a  homogeneous  circu¬ 
lar  conductor  cross  section,  phase  lag.  Problem  radius  R  = 
5/jm  and  co=  1GHz. 


Im(A.) 


Figure  7:  2D,  split-conductor  problem,  the  imaginary  part 
of  the  solution. 


As  expected,  the  total  current  density  is  maximum 
and  equal  to  the  dc  value  at  the  boundary,  and 
reaches  its  minimum  at  the  center  of  the  cross  sec¬ 
tion.  The  characteristic  skin-depth  decay  scale  is  well 
in  evidence  in  Fig.  (4).  In  addition,  the  expected 
maximum  phase  lag  occurs  at  p  =  0  in  Fig.  (5).  Table 
(1)  summarizes  the  results  for  the  skin-effect  prob¬ 
lem,  while  Table  (2)  shows  the  results  for  the  fre¬ 
quency-dependent  self  impedance  of  a  cross  section 
of  radius  1.0  fjm  at  frequencies  of  1  GHz,  5  GHz  and 
10  GHz.  As  expected,  both  the  frequency-dependent 
inductance  and  frequency  dependent  inductance  in¬ 
creases  with  frequency.  For  extracting  impedance,  a 
total  of  only  1,000  RWs  were  performed  per  points.  It 
can  be  seen  from  Table  (2),  that  the  error  in  the  esti¬ 
mate  of  frequency-dependent  resistance  and  inductive 
impedance  is  around  1  percent  in  all  three  cases.  Ta¬ 
ble  (3)  summarizes  the  results  for  the  heterogeneous 
problem,  while  Figures  (6)  and  (7)  illustrates  the  re¬ 
sults  for  the  same.  Again,  excellent  conformity  was 
obtained  between  the  analytical  and  RW  solutions. 


Re(Az) 


Figure  6:  2D,  split-conductor  problem,  the  real  part  of  the 
solution. 


Frequency 

(GHz) 

Random 

Walks 

Per  So¬ 
lution 
Point 

Mean  Abso¬ 
lute  Error 
For  Relative 
Magnitude 

Mean  Abso¬ 
lute  Error 
For  Relative 
Phase 

1 

20000 

0.001 

on  a  solution 
range 

(0-42 -1-0). 

0.012 

on  a  solution 

range 

(-1-9-0). 

Table  1:  Numerical  results  for  the  2D  skin-effect  problem 
in  a  circular  conducting  cross  section. 


Frequency 

(GHz) 

Time 

(Sec¬ 

onds) 

Ana¬ 

lytical 

Result 

(Q/m) 

RW 

Result 

(Q/m) 

Error 

(Q/m) 

1 

20 

5735 

5738 

3 

+ 

+ 

- 

314i 

312i 

2i 

5 

30 

5870 

5917 

47 

+ 

+ 

- 

1552i 

1534i 

18i 

10 

45 

6262 

6315 

53 

+ 

+ 

- 

2997i 

2962i 

35i 

Table  2:  Numerical  results  for  the  frequency-dependent 
self-impedance  of  a  conducting  circular  cross  section. 
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rithm  can  be  used  to  extract  electrical  parameters 
such  as  frequency-dependent  impedance.  We  believe 
that  with  additional  development,  the  algorithm  may 
prove  useful  for  electromagnetic  analysis  of  complex, 
multilevel  IC-interconnect  structures. 

Importantly,  the  algorithm  is  fully  parallel.  Thus,  we 
expect  significant  performance  acceleration  in  any 
future  parallel  software  or  hardware  implementation. 
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We  will  finish  our  discussion  in  this  section  by  mak¬ 
ing  a  few  comments  on  the  accuracy  of  the  solution, 
time  and  memory  requirement.  We  have  already  ob¬ 
served  the  close  agreement  of  the  RW  results  with 
that  of  known  analytical  solutions.  The  accuracy  of 
the  solution  of  the  solution  can  be  enhanced  by  sim¬ 
ply  increasing  the  number  of  RWs  as  the  error  is  pro¬ 
portional  to  1  /y[N~,  N  being  the  number  of  RWs. 
This  particular  fact  is  a  direct  consequence  of  Central 
Limit  Theorem  [14].  The  memory  requirements  for 
this  technique  are  low  as  this  approach  does  not  re¬ 
quire  any  numerical  meshing.  The  time  requirements 
of  this  algorithm  can  be  further  reduced  by  the  use  of 
variance-reduction  techniques  [2]  and  by  paralleliza¬ 
tion.  We  plan  to  investigate  all  these  issues  in  detail 
after  we  have  applied  our  algorithm  to  more  compli¬ 
cated  structures. 


VI.  CONCLUSION 

We  have  presented  the  theoretical  basis  of  a  novel 
floating  RW  algorithm  for  solving  the  2D  Maxwell- 
Helmholtz  equation.  The  algorithm  employs  iterative 
perturbation  theory.  We  have,  as  well,  verified  the 
algorithm’s  integrity  by  applying  it  to  a  homogeneous 
and  a  heterogeneous  problem,  possessing  analytical 
solutions.  The  applicability  to  heterogeneous  prob¬ 
lems  is  a  significant  improvement  on  existing  RW 
algorithms,  an  application  we  wish  to  explore  further 
in  our  future  work.  Our  algorithm  can  be  readily  ex¬ 
tended  to  multi-conductor  systems  in  full  3D.  In  this 
work,  we  have  further  demonstrated  that  this  algo- 
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Abstract:  It  is  usually  ignored  that 
the  transmission  line  parameters 
extracted  by  conventional  transmission 
line  equations  are  variable  with  the 
internal  impedance  of  excitation  source 
and  the  loading  of  the  line.  To  obtain 
the  correct  parameters,  the  match  for 
both  excitation  and  loading  has  to  be 
applied.  Problem  is  that  exact  matches 
for  some  extreme  cases  are  not  easy  to 
realize.  To  solve  this  problem,  in  this 
paper,  generalized  transmission-line 
equations  are  implemented  to  the 
transmission-line  parameter  extraction 
with  the  cooperation  of  full-wave  solver 
Zeland  IE3D.  The  parameters  extracted 
by  the  generalized  transmission  line 
equations  are  invariant  to  both  excitation 
and  loading.  Except  for  this,  the  local 
radiation  parameters  generated  from  the 
transmission  line  discontinuities  can 
also  be  found. 

Keywords:  generalized  transmission 
line  equation,  parameters  extraction,  and 
invariance. 


I.  INTRODUCTION 

For  a  simply  shaped  uniform 
transmission-line,  its  per-unit-length 
parameters  are  easily  found  by 
analytical  formula  and  the 
transmission-line  equations  can  be  used 


to  solve  the  line’s  problem  with  arbitrary 
excitation  and  load.  For  a  complicatedly 
shaped  non-uniform  transmission-line, 
however,  its  per-unit-length  parameters 
are  hard  to  be  determined  by  analytical 
formula,  and  have  to  be  extracted  by 
numerical  methods.  When  the 
conventional  transmission-line 

equations  are  used  to  extract  line’s 
per-unit-length  parameters,  what  we 
find  is  that  the  extracted  parameters  vary 
with  the  excitation  internal  impedance 
and  the  line  load.  The  reason  for  this  is 
that  the  derivation  of  the  conventional 
transmission-line  equations  is  based  on 
the  assumption  of  an  infinite  line  or 
non-reflection  that  is  equivalent  to  the 
matching  conditions  for  both  excitation 
and  load.  In  other  words,  the  extracted 
parameters  are  incorrect  if  the  matching 
condition  for  both  excitation  and  load  is 
not  satisfied.  It  should  be  emphasized 
that  the  matched  condition  has  been 
used  in  finding  correct  equation 
coefficients  only.  After  the  parameters 
(equation  coefficients)  are  determined, 
the  equations  can  be  used  to  solve  the 
line’s  problem  with  arbitrary  excitation 
and  load.  Unfortunately,  the  exact 
matches  for  some  extreme  cases,  such  as 
a  dipole  antenna,  are  hard  to  be  obtained. 
To  solve  this  problem,  in  this  paper, 
generalized  transmission-line  equations 
are  implemented  to  the  transmission-line 
parameter  extraction.  Since  the 
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derivation  of  the  generalized 
transmission-line  equations  is  based  on 
the  assumption  of  the  finite  rather  than 
infinite  line,  the  parameters  extracted  by 
the  generalized  transmission-line 
equations  are  invariance  to  both 
excitation  and  load.  In  comparison  with 
the  derivation  of  the  conventional 
transmission-line  equations,  the 
generalized  transmission-line  equations 
have  introduced  two  new  terms,  the 
dependent  voltage  source  aV  and 
dependent  current  source  pi ,  which  are 
interpreted  as  the  local  radiation  sources. 
Here,  V  and  /  stand  for  the  voltage  and 
current  at  each  discrete  segment  of  the 
transmission-line,  while  a  and  /?  are 
the  coefficients  of  V  and  I  , 
respectively.  In  order  to  extract  the  line 
parameters  by  using  the  generalized 
transmission  line  equations,  two  pairs  of 
voltage  and  current  solutions  along  the 
line  have  to  be  used.  By  means  of  the 
MoM  commercial  software  like  Zeland 
IE3D,  Sonnet,  and  Ensemble,  the  two 
pairs  of  solutions  for  two  different  loads 
can  be  obtained.  Thus,  the  parameters 
can  be  found  by  substituting  these  two 
pairs  of  solutions  into  the  generalized 
transmission  line  equations.  The 
parameters  extracted  by  the  generalized 
equations  are  invariant  to  both  internal 
source  impedance  and  loads.  To  show 
the  invariance  of  the  extracted 
parameters  to  the  excitation  and  load,  in 
this  paper,  two  numerical  examples  of 
microstrip  line  structures  are  given. 
Except  for  this,  the  local  radiation 
effects  from  the  discontinuities  can  also 
be  obtained  when  the  generalized 
equations  are  applied  to  non-uniform 
transmission  line  structures. 


II.  IMPLEMENTATION  OF 
GENERALIZED  EQUATIONS 
INTO  PARAMETER  EXTRACTION 


As  is  well  known,  the  following 
conventional  transmission  line  equations 
have  been  derived  from  Kirchhoff’s 
laws  on  a  basis  of  an  infinite-length  line 
and  TEM  mode  assumption  [2] 


dV(z) 

dz 

dl(z) 

dz 


■  =  -Z(z)-/(z) 
=  -T(z)-V(z) 


0) 


where  Z  =  jo)L+R  and 

Y  =  jcoC  +  G  are  per-unit-length  series 
impedance  and  shunt  admittance, 
respectively.  To  let  the  equations  be  used, 
the  per-unit-length  line  parameters 
L  ,R ,  C  and  G  (i.e.,  coefficients  of 
the  equations)  must  be  found  by  using 
either  analytical  formula  or  numerically 
extracted  technique.  In  the  past  years, 
almost  all  literatures  resorted  to  directly 
solving  Maxwell’s  equations  to  find  the 
parameters  under  static  and  quasi-static 
field  assumption  [2]  -  [5].  However, 
when  equations  (1)  are  used  in  the 
parameter  extraction,  one  of  things, 
which  is  easily  ignored,  is  that  the 
extracted  parameters  are  directly 
dependant  on  internal  impedance  and 
line’s  load  impedance.  In  other  words,  to 
obtain  correct  parameters,  both 
excitation  and  loading  matches  have  to 
be  imposed  so  that  the  assumption  of  the 
infinite  line  can  be  satisfied.  The 
problem  is  that  exact  matches  at  both 
ends  of  the  line  are  not  easy  to  realize 
because  the  characteristic  impedance  of 
the  line  is  unknown.  To  solve  this 
problem,  the  following  generalized 
transmission  line  equations  derived  from 
a  finite  line  [6] 
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—  =  -Z(/)/(/)  +  or(/)V(/) 

dl  (2) 

^-  =  -Y(l)V(l)+/3(l)I(l) 

al 

are  implemented.  For  a  lossless  line,  the 
distributed  parameters  Z(l )  and  Y(l) 
can  be  written  as  jcoL(J)  and  j(oC{l) , 
and  the  corresponding  definitions  are 
exactly  like  those  in  the  conventional 
transmission-line  equations. 

Comparing  the  conventional 
transmission-line  equations  (1),  the 
generalized  equations  introduce  two 
new  parameters  ail)  and  /?(/) ,  which 
can  be  interpreted  as  the  coefficients  of 
dependent  voltage  source  and  dependent 
current  source  in  circuit  theory,  or  be 
interpreted  as  the  coefficients  of  local 
radiations  between  discontinuities  of 
neighboring  segments  in  field  theory. 
When  the  transmission-line  is  a  uniform 
structure,  the  values  of  ail)  and  /?(/)  are 
to  be  zero  so  that  the  generalized 
equations  have  the  same  formulations  as 
the  conventional  uniform  line  equations. 
At  first  glance,  the  generalized 
equations  seem  to  be  more  complex  than 
the  conventional  equations  because 
there  are  four  parameters  to  be  extracted 
at  each  discrete  segment.  In  fact,  as 
stated  above,  the  two  additional 
parameters  make  the  extracted 
parameters  be  invariant  to  line’s 
excitation  and  load  so  that  the  correct 
line’s  parameters  can  easily  extracted. 
Since  there  are  four  line  parameters 
Z(/) ,  Y(l) ,  ail)  and  /?(/)  at  each 
discrete  segment  in  generalized 
equations  (2),  we  need  to  use  two  pairs 
of  the  distributed  voltage  and  current 
solutions  along  the  transmission-line  to 
determine  the  line  four  parameters.  Now 
most  full-wave  solver  tools  like  Zeland 


IE3D,  Sonnet,  and  Ensemble  can  be 
used  for  this  purpose. 

m.  NUMERICAL  EXAMPLES 

To  show  how  to  implement  the 
generalized  equations  to  the  line 
parameter  extraction,  the  following  two 
examples  are  presented.  One  is  a 
finite-length  uniform  microstrip  line. 
Another  is  a  microstrip  bend  with  two 
arms.  For  these  two  structures,  the 
relative  dielectric  constant  er  of  the 
substrate  is  9.8,  the  height  h  between  the 
metal  strip  and  ground  is  0.635mm,  the 
thickness  t  of  the  metal  strip  is  2/jm, 
and  the  width  w  of  the  metal  strip  is 
0.6mm.  For  simplicity,  the  metal  are 
supposed  to  be  lossless.  The  full-wave 
simulation  tool,  Zeland  IE3D  Software 
[7],  is  used  to  compute  the  distributed 
voltages  and  currents  along  the 
microstrip  structures. 

For  the  first  example,  the  line 
parameters  UJ),  C(J),  ail)  and 
/?(/)  for  three  cases  are  extracted  by  the 
generalized  transmission  line  equations. 
The  first  case  is  for  short  and  open  loads 
(Al).  The  second  case  is  for  the  loads 
with  complex  numbers  of  20+20j  and 
lOO+lOOj  (Bl).  The  third  case  is  for  the 
complex  number  internal  impedance  of 
20+20j  and  lOO+lOOj  for  excitation 
source  (Cl).  The  internal  source 
impedances  for  the  first  two  cases  are 
50 Cl,  and  terminated  load  for  the  last 
case  is  50  Q .  All  of  the  parameters  Ul), 
C(l),  a(l)  and  J3il)  along  the 
uniform  microstrip  line  found  by  the 
equations  (2)  are  shown  in  Figures  1-4, 
respectively.  It  can  be  seen  that  the 
extracted  parameters  are  almost  the 
same  for  the  three  cases.  This  implies 
that  the  line  parameters  extracted  by  the 
generalized  equations  are  invariant  to 
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Fig.  1  Comparison  of  inductance  obtained  from 
generalized  equations  and  traditional 
equations  at  1GHz. 


l(mm) 

Fig.  3  Distributions  of  alpha  obtained  from 
generalized  equations  at  1GHz. 

both  internal  source  impedance  and 
termination  loads.  While  the  same  cases 
are  applied  to  the  conventional 
equations  (1)  (corresponding  to  cases 
A2,  B2  and  C2),  the  extracted  L(/)  and 
C(/)  are  variant  with  both  excitation  and 
load  except  for  the  excitation  and  load 
matching  conditions  (case  F),  as  show  in 
Figs.  1  and  2.  In  other  words,  the  values 
of  L(l)  and  C(/)  extracted  by  the 
conventional  equations  (1)  at  matching 


-30.0  -20.0  -10.0  0.0  10.0  20.0  30.0 

l(mm) 


Fig.  2  Comparison  of  capacitance  obtained  from 
generalized  equations  and  traditional 
equations  at  1GHz. 


Fig.  4  Distributions  of  beta  obtained  from 
generalized  equations  at  1GHz. 


conditions  (case  F)  are  of  a  good 
agreement  with  those  extracted  by 
generalized  equations  (2)  and  calculated 
by  analytical  formulae  [3],  but  the 
values  of  L{t)  and  C(/)  extracted  by  the 
conventional  equations  (1)  at 
mismatching  conditions  are  at  a  great 
difference  from  those  extracted  with 
generalized  equations  (2).  In  addition, 
we  find  that  the  local  radiation 
coefficients  #(/)  and  /?(/) ,  as  shown 
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in  Fig.  3  and  Fig.  4,  are  near  to  zero 
except  for  in  the  vicinity  of  both  line 
ends.  The  non-zero  #(/)  and  /?(/)  in 
the  vicinity  of  both  line  ends  are 
reasonable  because  the  line  ends  for  a 
finite-length  transmission  line  itself  just 
are  discontinuous  places  of  the  line.  For 
second  example,  the  microstrip  bend  is 
not  uniform  structure,  so  the 


0.60  - 


l(mm) 

Fig.  5  Comparison  of  inductance  obtained 
from  generalized  equations  with 
different  loads  at  1GHz. 


40.0 - 
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Fig.  7  Comparison  of  alpha  obtained 
from  generalized  equations  with 
different  loads  at  1GHz. 


conventional  transmission  line  equations 
could  not  be  directly  applied  to  it. 
However,  the  generalized  equations  can 
be  directly  used.  Two  cases  are 
calculated  for  this  bend  structure.  Both 
cases  have  the  excitation  source  of  50  Q . 
The  first  case  is  for  short  and  open  loads 
(Al).  The  second  case  is  for  20+20j  and 
100+lOOj  loads  (Bl). 


165.0  .- 
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l(mm) 

Fig.  6  Comparison  of  capacitance  obtained  from 
generalized  equations  with  different  loads 
at  1GHz. 


0.20 


l(mrri) 


Fig.  8  Comparison  of  beta  obtained 
from  generalized  equations  with 
different  loads  at  1GHz. 
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The  line  parameters  L(l),  C(l),  a{l) 
and  /?(/)  are  extracted  by  the 
generalized  equations  (2).  Figures  5-8 
show  the  results  of  UJ),  C(I),  a{l)  and 
/?(/)  along  the  whole  bend  structure. 
For  the  second  case,  the  line 
characteristic  impedance  is  far  away 
from  that  of  traveling  waves.  However, 
like  the  first  example,  the  two  cases  of 
parameters  are  also  invariant  to  the  port 
conditions.  In  addition,  the 
characteristics  of  the  bend  discontinuity 
can  be  also  obtained,  as  shown  in  Fig.  7 
and  Fig.  8,  where  the  larger  values  of 
a(I )  and  /?(/)  appear  around  the 
comer  of  bend  while  near-zero  values 
occur  at  the  flat  part  of  two  arms. 

IV.  CONCLUSION 

In  this  paper,  the  generalized 
transmission-line  equations  are 
implemented  to  extract  the  transmission 
line  parameters  by  means  of  numerical 
simulation  tools  like  Zeland  IE3D 
Software.  The  distinguished  property  for 
the  general  equations  is  that  the 
extracted  line  parameters  are  invariant  to 
both  internal  source  impedance  and 
terminated  loads.  Besides,  the  local 
radiation  effects  from  line 
discontinuities  can  also  be  obtained 
during  the  process  of  the  parameter 
extraction. 
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